trsing’s diary

勉強、読んだ本、仕事で調べたこととかのメモ

PRML 4.3.4~4.3.6

4.3.4 多クラスロジスティック回帰

多クラスの分布に対して事後確率p(C_{k}|\phi)は特徴変数の線形関数のソフトマックス変換で与えられる。この節では最尤法を用いて直接パラメータ\{w_{k}\}を決定する。

パラメータベクトルwに関する線形回帰モデルの対数尤度関数の微分は誤差y_{n}-t_{n}と特徴ベクトル\phi_{n}との積になる。これは今回に限らずもっと一般的な結果。(4.3.6節で見る)

多クラス問題用のIRLSアルゴリズムを得るために必要なヘッセ行列のブロックj,kは(4.110)となる。

(4.106):
f:id:trsing:20181213223740j:plain

(4.109): f:id:trsing:20181213225102j:plain

(4.110):
f:id:trsing:20181213225930j:plain

4.3.5 プロビット回帰

指数型分布族以外のタイプの識別確率モデルについて調べる。 ただし、2クラスの場合で(4.111)の範囲で議論を行う。

プロビット回帰について。プロビット活性化関数に基づく一般化線形モデル。S字形。ロジスティック回帰より外れ値に対して敏感。

(4.116):
f:id:trsing:20181213233256j:plain

4.3.6 正準連結関数

これまでにいくつかのモデルで パラメータベクトルwに関する誤差関数の微分が「誤差」y_{n}-t_{n}と特徴ベクトル\phi_{n}との積になるのを見てきた。 これは正準連結関数として知られている関数を活性化関数に選び、指数型分布族の中から目的変数に対する条件付確率分布を選択することから得られる一般的な結果であることを示す。

(4.119) :
f:id:trsing:20181213233350j:plain

(4.122):

\frac{d \eta_{n}}{d y_{n}}=\frac{d \varphi(y_{n})}{d y_{n}}=\varphi'(y_{n})\\
\frac{d y_{n}}{d a_{n}}=\frac{d f(w^{T}\phi_{n})}{d a_{n}}=\frac{d f(a_{n})}{d a_{n}}=f'(a_{n})

(4.123): 
f^{-1}(y)=\varphi(y)\longrightarrow y=f(\varphi(y))\longrightarrow f'(\varphi(y))\varphi'(y)=\frac{d f(\varphi(y))}{d \varphi(y)}\frac{d\varphi(y)}{dy}=\frac{dy}{d\varphi(y)}\frac{d\varphi(y)}{dy}=1\\
a=f^{-1}(y)=\varphi(y)\\
f'(a)\varphi'(y)=\frac{d f(a)}{da}\frac{df(\varphi(y))}{dy}=\frac{d f(\varphi(y))}{d \varphi(y)}\frac{d\varphi(y)}{dy}=1

PRML 4.3~4.3.3

4.3 確率的識別モデル

一般化線形モデルの関数形式を陽に仮定し、パラメータを直接決定する方法。
条件付確率密度分布p(C_{k}|x)を通じて定義される尤度関数を最大化する(識別学習の一形態)。

利点
 決めるべき適応パラメータが少ない。
 真の確率分布をうまく近似できない場合でもよい性能を示す場合がある。

4.3.1 固定基底関数

基底関数ベクトル\phi(x)を使って入力の非線形変換を行っておけばこれまでのアルゴリズムを同じように適用できる。

決定境界は特徴空間\phi(x)において線形だが、もとの入力空間xにおいては非線形の決定境界となる。特徴空間\phi(x)で線形分離可能であるクラスが、もとの観測空間xで線形分離可能である必要はない。

線形変換\phi(x)はクラス間の重なりを取り去ることはできないが、非線形性を適切に選択すれば、事後確率のモデル化が簡単になる

4.3.2 ロジスティック回帰

(4.87):ロジスティック回帰。分類のためのモデル。
M次元特徴空間\phiにおいて調整可能なパラメータ数はM個。 ガウス分布の場合次元Mの二乗で増加するため、ロジスティック回帰のほうが有利。

パラメータを最尤法を用いて決定するために交差エントロピー誤差関数を与え、これを微分すると(4.91)。これは(3.13)と同じ形。ただしロジスティックシグモイド関数非線形性のため最尤解を解析的に導出できない。

線形分離可能なデータ集合に対しては過学習を起こすので注意。

ガウス分布のパラメータ数 : 各クラスの平均(\mu_{1},\mu_{2})に対してM個ずつ。共分散行列は共通、対象なのでM(M+1)/2

(4.88):
(4.61)を使用して

\frac{d\sigma}{da}=\frac{\exp(-a)}{(1+\exp(-a))^{2}}=\frac{\frac{1-\sigma}{\sigma}}{(1+\frac{1-\sigma}{\sigma})^{2}}=\sigma(1-\sigma)

(4.91):
(4.88)を使用して

\frac{\partial}{\partial w}\ln\sigma(w^{T}\phi_{n})=
\frac{\partial}{\partial \sigma}\ln\sigma
\frac{\partial \sigma}{\partial a}
\frac{\partial a}{\partial w}
=\frac{1}{\sigma}\sigma(1-\sigma)\phi_{n}
=(1-\sigma)\phi_{n}

\frac{\partial}{\partial w}\ln(1-\sigma(w^{T}\phi_{n}))=
\frac{\partial}{\partial 1-\sigma}\ln(1-\sigma)
\frac{\partial (1-\sigma)}{\partial a}
\frac{\partial a}{\partial w}
=\frac{1}{1-\sigma}(-\sigma)(1-\sigma)\phi_{n}
=-\sigma\phi_{n}

4.3.3 反復最重み付け最小二乗

ロジスティック回帰はロジスティックシグモイド関数非線形性により、最尤解を解析的に導出することはできない。
しかし、誤差関数は凸関数なので、唯一の最小解をもち、ニュートン‐ラフソン法に基づく反復最適化手順を用いて最小化できる。

線形回帰モデル(3.3)にニュートン‐ラフソン法を適用すると(4.95)となりこれは標準的な最小二乗解となる。

ロジスティック回帰における交差エントロピー誤差関数にニュートン‐ラフソン法を適用すると(4.99)となる。パラメータベクトルwが更新されるたびに重みづけ行列Rの再計算が必要なため、反復最重み付け最小二乗法として知られている。

(4.101):

t\in\{0,1\}\\
p(t|w)=y^{t}(1-y)^{1-t}\\
E[t]=\int tp(t|w) dt=y

(4.103):

\frac{da}{d y}=\frac{da}{d \sigma}=\frac{1}{\sigma(1-\sigma)}=\frac{1}{y(1-y)}

PRML 4.1.7~4.2.4

4.1.7 パーセプトロンアルゴリズム

誤分類された場合w^{T}\phi_{n}t_{n}は負になる。
判定基準が正負なので定数を掛けても変化しない

(\phi_{n}t_{n})^{T}\phi_{n}t_{n}=||\phi_{n}t_{n}||^{2}

4.2 確率的生成モデル

線形決定境界をどのように生成するか示す。
ロジスティックシグモイド関数、ソフトマックス関数の導入。ソフトマックス関数はロジスティックシグモイド関数の多クラスへの一般化とみなすことができる。

(4.61)

\exp(-a)=\frac{1-\sigma(a)}{\sigma(a)}\\
a=\ln\frac{\sigma(a)}{1-\sigma(a)}

4.2.1 連続値入力

仮定
・クラスの条件付確率密度がガウス分布であると仮定する。
・すべてのクラスが同じ共分散行列を共有すると仮定する。
を置くと
・2クラスの場合、決定境界は入力関数xで線形になる。
・Kクラス分類の一般的な場合もxの線形関数として定義される。
共分散行列が共通ではない場合、線形決定境界は2次となる

(4.65)

p(C_{1}|x)=\sigma(a)\\
a=\ln\frac{p(x|C_{1})p(C_{1})}{p(x|C_{2})p(C_{2})}
p(x|C_{1})p(x|C_{2})の共通項( \frac{1}{(2\pi)^{D/2}}\frac{1}{|\Sigma|^{1/2}}
\exp(-\frac{1}{2}x^{T}\Sigma^{-1}x)) をキャンセル。

(4.68)
(4.62)より分子分母の共通項(\frac{1}{(2\pi)^{D/2}}\frac{1}{|\Sigma|^{1/2}}
,\,\exp(-\frac{1}{2}x^{T}\Sigma^{-1}x))をキャンセルする。

4.2.2 最尤解

データ集合\{x_{n},t_{n}\}から事前分布p(C_{k})と各クラスの条件付確率密度(p(x_{n}|C_{k})ガウス分布)のパラメータを求める。

まぎらわしいので\piじゃなくて\pi'と書く。
(4.71)の対数尤度関数。
\ln p(t,X|\pi',\mu_{1},\mu_{2},\Sigma)=\sum[t_{n}\ln\pi'+(1-t_{n})\ln(1-\pi')-\frac{D}{2}\ln2\pi-\frac{1}{2}\ln|\Sigma|\\
\,\,-\frac{t_{n}}{2}(x_{n}-\mu_{1})^{T}\Sigma^{-1}(x_{n}-\mu_{1})-\frac{1-t_{n}}{2}(x_{n}-\mu_{2})^{T}\Sigma^{-1}(x_{n}-\mu_{2})]

\pi'微分
\sum(\frac{t_{n}}{\pi'}-\frac{1-t_{n}}{1-\pi'})=\sum\frac{t_{n}-\pi'}{\pi'(1-\pi')}=0\\
\sum\pi'=\sum t_{n}\\
\pi'=\frac{1}{N}\sum t_{n}

\mu_{1}微分

\sum t_{n}\Sigma^{-1}(x_{n}-\mu_{1})=0\\
\sum t_{n}\mu_{1}=\sum t_{n}x_{n}\\
\mu_{1}=\frac{1}{N_{1}}\sum t_{n} x_{n}

\mu_{2}微分

\sum (1-t_{n}) \Sigma^{-1}(x_{n}-\mu_{2})=0\\
\sum (1-t_{n}) \mu_{2}=\sum (1-t_{n})x_{n}\\
\mu_{2}=\frac{1}{N_{2}}\sum (1-t_{n}) x_{n}

\Sigmaについて整理すると

\sum t_{n}(x_{n}-\mu_{1})^{T}\Sigma^{-1}(x_{n}-\mu_{1})\\
\sum_{n \in C_{1}} Tr[(x_{n}-\mu_{1})^{T}\Sigma^{-1}(x_{n}-\mu_{1})]\\
\sum_{n \in C_{1}} Tr[\Sigma^{-1}(x_{n}-\mu_{1})(x_{n}-\mu_{1})^{T}]\\
Tr[\Sigma^{-1}\sum_{n \in C_{1}}(x_{n}-\mu_{1})(x_{n}-\mu_{1})^{T}]
より(4.77)

\Sigma微分すると

\frac{\partial}{\partial X} \ln|X|=(X^{-1})^{T}\\
\frac{\partial}{\partial X} Tr[AX^{-1}B]=-(X^{-1}BAX^{-1})^{T}
より

-N/2(\Sigma^{-1})^{T}+N/2(\Sigma^{-1}S\Sigma{-1})^{T}=0\\
(\Sigma^{-1}S)^{T}=I\\
\Sigma=S

4.2.3 離散特徴

特徴が離散値x_{i}の場合。

4.2.4 指数型分布

これまでの話( ガウス分布と離散値入力の場合、クラスの事後確率はロジスティックシグモイド関数またはソフトマックス活性化関数の一般化線形モデルで与えられる) は仮定(クラスの条件付確率密度p(x|C_{k})が指数型分布族のメンバーである)から得られる結果の特殊な場合。

指数型分布族(4.83)にてu(x)=xとなるような分布の場合、クラスの事後確率はxの線形関数のロジスティックシグモイド関数、ソフトマックス関数で与えられる。

(4.86)
(4.62)から共通項(\frac{1}{s}h(\frac{1}{s}x))を消してる

PRML 4.1.3~4.1.6

4.1.3 分類における最小二乗

(4.15)を微分

\frac{\partial}{\partial X}Tr[(AXB+C)(AXB+C)^{T}]
=2A^{T}(AXB+C)B^{T},\\
Tr[(AXB+C)(AXB+C)^{T}]=Tr[(AXB+C)^{T}(AXB+C)]
より
\frac{\partial}{\partial W}E_{D}(W)=2X^{T}(XW-T)

4.1.4 フィッシャーの線形判別

(演習4.4) w^{T}(m_{2}-m_{1})+\lambda(w^{T}w-1)w微分するとm_{2}-m_{1}+\lambda 2w

(4.25)から(4.26)

(y_{n}-m_{k})^{2}=(w^{T}x_{n}-w^{T}m_{k})(w^{T}x_{n}-w^{T}m_{k})^{T}=w^{T}(x_{n}-m_{k})(x_{n}-m_{k})^{T}w
(m_{2}-m_{1})^{2}も同様

(4.26)を微分

\frac{\partial}{\partial x}\frac{(Ax)^{T}Ax}{(Bx)^{T}Bx}=2\frac{A^{T}Ax}{(Bx)^{T}Bx}-2\frac{(x^{T}A^{T}Ax)B^{T}Bx}{(x^{T}B^{T}Bx)^{2}}
より

\frac{\partial}{\partial w}J(w)=2\frac{S_{B}w}{w^{T}S_{w}w}-2\frac{(w^{T}S_{B}w)S_{w}w}{(w^{T}S_{w}w)^{2}}

S_{B}w=(m_{2}-m_{1})(m_{2}-m_{1})^{T}w(m_{2}-m_{1})^{T}wスカラーなのでS_{B}w(m_{2}-m_{1})と同じ方向を持つ。

4.1.5 最小二乗との関係

$$ \sum_{n=1}^{N}x_{n} = \sum_{k=1}^{K} \sum_{n \in C_{k}} x_{n}=\sum_{k=1}^{K} N_{k}m_{k} $$ に注意すると
(4.32)=w^{T}\sum x_{n}+Nw_{0}-\sum t_{n}=N_{1}m_{1}+N_{2}m_{2}+N w_{0}
(4.33)=\sum(w^{T}x_{n})x_{n}-(w^{T}m)\sum x_{n}-\sum t_{n}x_{n}\\
=\sum (x_{n}^{T}w)x_{n}-(m^{T}w)(N_{1}m_{1}+N_{2}m_{2}) -N/N_{1}\sum x_{n}+N/N_{2}\sum x_{n}\\
=\sum x_{n}(x_{n}^{T}w)-(N_{1}m_{1}+N_{2}m_{2})(m^{T}w)-N(1/N_{1}\sum x_{n}-1/N_{2}\sum x_{n})\\
=(\sum x_{n}x_{n}^{T})w-(N_{1}m_{1}+N_{2}m_{2})1/N(N_{1}m_{1}+N_{2}m_{2})^{T}w-N(m_{1}-m_{2})\\
=0
(S_{W}+N_{1}N_{2}/NS_{B})=\sum
(x_{n}x_{n}^{T}-m_{1}x_{n}^{T}-x_{n}m_{1}^{T}+m_{1}m_{1}^{T})\\
\hspace{100pt}+\sum
(x_{n}x_{n}^{T}-m_{2}x_{n}^{T}-x_{n}m_{2}^{T}+m_{2}m_{2}^{T})\\
\hspace{100pt}+N_{1}N_{2}/N(m_{2}m_{2}^{T}-m_{1}m_{2}^{T}-m_{2}m_{1}^{T}+m_{1}m_{1}^{T})\\
\hspace{85pt}=(\sum
x_{n}x_{n}^{T})-m_{1}N_{1}m_{1}^{T}-N_{1}m_{1}m_{1}^{T}+N_{1}m_{1}m_{1}^{T}\\
\hspace{100pt}+(\sum
x_{n}x_{n}^{T})-m_{2}N_{2}m_{2}^{T}-N_{2}m_{2}m_{2}^{T}+N_{2}m_{2}m_{2}^{T}\\
\hspace{100pt}+N_{1}N_{2}/N(m_{2}m_{2}^{T}-m_{1}m_{2}^{T}-m_{2}m_{1}^{T}+m_{1}m_{1}^{T})\\
\hspace{85pt}=(\sum x_{n}x_{n}^{T})-N_{1}m_{1}m_{1}^{T}-N_{2}m_{2}m_{2}^{T}\\
\hspace{100pt}+N_{1}N_{2}/N(m_{2}m_{2}^{T}-m_{1}m_{2}^{T}-m_{2}m_{1}^{T}+m_{1}m_{1}^{T})\\
\hspace{85pt}=(\sum
x_{n}x_{n}^{T})
-1/N
(N_{1}(N_{1}+N_{2})m_{1}m_{1}^{T}+N_{2}(N_{1}+N_{2})m_{2}m_{2}^{T}
-N_{1}N_{2}(m_{2}m_{2}^{T}-m_{1}m_{2}^{T}-m_{2}m_{1}^{T}+m_{1}m_{1}^{T}))\\
\hspace{85pt}=(\sum
x_{n}x_{n}^{T})
-1/N
(N_{1}^{2}m_{1}m_{1}^{T}+N_{1}N_{2}(m_{1}m_{2}^{T}+m_{2}m_{1}^{T})+N_{2}^{2}m_{2}m_{2}^{T})\\
\hspace{85pt}=(\sum
x_{n}x_{n}^{T})
-1/N
(N_{1}m_{1}+N_{2}m_{2})(N_{1}m_{1}+N_{2}m_{2})^{T}\\
これより(4.33)と(4.37)は等しい。

4.1.6 多クラスにおけるフィッシャーの判別

$$ S_{T}=\sum_{n=1}^{N}(x_{n}-m)(x_{n}-m)^{T}=\sum_{k=1}^{K}\sum_{n \in C_{k}} (x_{n}-m)(x_{n}-m)^{T}\\= \sum_{k=1}^{K}\sum_{n \in C_{k}}( (x_{n}-m_{k})+(m_{k}-m) )( (x_{n}-m_{k})+(m_{k}-m) )^{T}\\= \sum_{k=1}^{K} \sum_{n \in C_{k}} [ (x_{n}-m_{k})(x_{n}-m_{k})^{T}+ (x_{n}-m_{k})(m_{k}-m)^{T}+ (m_{k}-m)(x_{n}-m_{k})^{T}+ (m_{k}-m)(m_{k}-m)^{T} ]\\= \sum_{k=1}^{K} [ S_{k}+N_{k}(m_{k}-m)(m_{k}-m)^{T}+ (\sum (x_{n}-m_{k}) )(m_{k}-m)^{T}+ (m_{k}-m)(\sum (x_{n}-m_{k}) )^{T} ]\\= S_{W}+S_{B}+ \sum_{k=1}^{K} [ (N_{k}m_{k}-N_{k}m_{k})(m_{k}-m)^{T}+ (m_{k}-m)(N_{k}m_{k}-N_{k}m_{k})^{T} ]\\= S_{W}+S_{B} $$

4~4.1.1

第四章 線形識別モデル

分類問題に対する3つのアプローチ

  • 各入力ベクトルxから直接クラスを推定する識別関数を構築する方法

  • 条件付き確率分布p(C_{k}|x)を直接モデル化する方法

  • クラスに対する事前確立p(C_{k})とともにp(x|C_{k})で与えられるクラスで条件づけされた確率密度を生成し、ベイズ定理を使って事後確率p(C_{k}|x)を計算する方法

4.1.1 2クラス

(4.5):xからwへの射影は||x||cos\thetaw^{T}x=||w||||x||cos\thetaより

(4.7):y(x)=w^{T}x_{\perp}+\frac{rw^{T}w}{||w||}=0+\frac{r||w||^{2}}{||w||}

PRMLメモ 3章 3.4~3.5.3

3.4 ベイズモデル比較

p(\mathcal{D}|\mathcal{M}_{i}):モデルエビデンス、または周辺尤度
p(\mathcal{D}|\mathcal{M}_{i})/p(\mathcal{D}|\mathcal{M}_{j})ベイズ因子

単純なモデルで得られるデータ集合\mathcal{D}は多様性に乏しい:例えばy=w_{0}であれば得られるデータ集合はw_{0} (+ランダムノイズ)。y=w_{0}+w_{1}xであれば得られるデータ集合は直線に乗るもの(+ランダムノイズ)。
複雑なモデルでは多様なデータ集合が得られるがそのため特定のデータ集合\mathcal{D}_{0}が得られる確率は低くなる。

(3.73)は対数を取ったベイズ因子の(モデル\mathcal{M}_{1}に従う)データ集合に関する期待値。E_{D}[ln \frac{p(\mathcal{D}|\mathcal{M}_{1})}{p(\mathcal{D}|\mathcal{M}_{2})}|\mathcal{M}_{1}]。これが>0より平均的にp(\mathcal{D}|\mathcal{M}_{1})>p(\mathcal{D}|\mathcal{M}_{2})

3.5 エビデンス近似

エビデンス近似:パラメータwだけに関して積分して周辺尤度関数を得る。これを最大にするように超パラメータの値を決める。
事前分布p(\alpha,\beta)が平坦→定数として扱う→p(t|\alpha,\beta)を最大にする\alpha,\betap(\alpha,\beta|t)の最大値が得られる。

w,\alpha,\betaに関して周辺化して予測分布を得る(3.74)
\alpha,\betaを固定してwを周辺化して予測分布を得る(3.75)
・事後分布p(\alpha,\beta|t)を最大にする\alpha,\betaを得たい。→(3.76)からp(t|\alpha,\beta)を最大にすることで得られる。(p(\alpha,\beta)は比較的平坦とする)

\alpha/\beta正則化パラメータと同様の働きをする:(3.55)参照

3.5.1 エビデンス関数の評価

エビデンス関数(3.77)の対数表現は(3.86)。導出はおおむね演習問題で。
(3.83):\frac{\partial^{2} E(w)}{\partial w_{i} \partial w_{j}}=\delta_{ij}\alpha+\beta\sum\phi_{i}(x_{n})\phi_{j}(x_{n})より
(3.86)は(3.72)に対応(\mathcal{D}tが対応?)

3.5.2 エビデンス関数の最大化

\alpha\sum\frac{1}{\lambda_{i}+\alpha}=\sum\frac{\lambda_{i}+\alpha-\lambda_{i}}{\lambda_{i}+\alpha}=M-\sum\frac{\lambda_{i}}{\lambda_{i}+\alpha}
\Phi^{T}\Phiが変化しないので\Phi^{T}\Phi固有ベクトルu_{i}も変化しない。
\frac{d}{d\beta}\lambda_{i}u_{i}=\frac{d}{d\beta}\beta\Phi^{T}\Phi u_{i}=\Phi^{T}\Phi u_{i}=\frac{\lambda_{i}}{\beta}u_{i}より\frac{d \lambda_{i}}{d\beta}=\frac{\lambda_{i}}{\beta}

3.5.3 有効パラメータ数

w_{ML}=(\Phi^{T}\Phi)^{-1}\Phi^{T}t
w_{MAP}=m_{N}=\beta(\alpha I+\beta\Phi^{T}\Phi)^{-1}\Phi^{T}t
\alpha=0ならw_{MAP}=w_{ML}\alphaが無限なら0。w_{MAP}0~w_{ML}間にある。
\alpha\lambda_{i}と比較して十分大きいならw_{MAP}i要素は0に近くなる。

(3.21):1/\beta_{ML}=1/N\sum\{t_{n}-w^{T}_{ML}\phi(x_{n})\}^{2}
(3.95):1/\beta=1/(N-\gamma)\sum\{t_{n}-m^{T}_{N}\phi(x_{n})\}^{2}

図3.16:交点が\gamma=\alpha m_{N}^Tm_{N}となる点。

(3.92)\alpha=\frac{\gamma}{m_N^{T}m_{N}}=\frac{\gamma}{2E_{W}(m_{N})}
(3.95)\beta=\frac{N-\gamma}{\sum\{t_{n}-m_{N}^{T}\phi(x_{n})\}^{2}}=\frac{N-\gamma}{2E_{D}(m_{N})}
N\gg Mの時(3.98),(3.99)

PRMLメモ 3章 図3.8、図3.9

図3.8と図3.9の図を描画するpythonのコード
参考(というかコメント以外ほぼ丸写し…)
qiita.com
コード見ると理解しやすいですね。

#https://qiita.com/naoya_t/items/80ea108cebc694f5cd63参照
from pylab import *

S=0.1
ALPHA=0.1
BETA=9

#(xs,ts):訓練データ。観測値と目標値。
#訓練データから事後分布の平均(3.53)と分散の逆数(3.54)を求める
#これらから予測分布(3.58)を求める
#予測分布(3.58)の平均と分散から図3.8を描画
#事後分布(3.49)からモデルパラメータwを決定して図3.9を描画
def sub(xs,ts):
    #s:標準偏差、mu:平均
    #\phi_j(x)作成する関数を返す。ガウス基底関数。
    def gaussian_basis_func(s,mu):
        return lambda x:exp(-(x - mu)**2/(2 * s**2))
    
    #s:標準偏差、xs:サイズ
    #基底関数\phi(x)を作成する関数を返す。
    def gaussian_basis_funcs(s,xs):
        return [gaussian_basis_func(s,mu) for mu in xs]
    
    xs9=arange(0,1.01,0.125)#9個
    bases=gaussian_basis_funcs(S,xs9)
    
    N=size(xs)
    M=size(bases)
    #\phi(x)作成
    def Phi(x):
        return array([basis(x) for basis in bases])
    
    PHI=array(list(map(Phi,xs)))#\Phiを作成
    PHI.resize(N,M)

    #重みの事後分布(3.49)の平均、分散、予測分布(3.58)の平均と分散を返す関数を返す
    def predictive_dist_func(alpha, beta):#wの平均m_N、分散S_N^{-1}、推定値の関数を返す
        #(3.49)の平均(m_N)と分散の逆行列(S_N^{-1})を作成
        S_N_inv=alpha*eye(M)+beta*dot(PHI.T,PHI)#(3.54):dotは積。.Tは転値
        m_N=beta*solve(S_N_inv,dot(PHI.T,ts))#(3.53):solve(A,b)はAx=bの解(A^(-1)b)
    
        #入力xに対応する予測分布(3.58)の平均と分散を返す
        def func(x):
            Phi_x=Phi(x)#\phi(x)
            mu=dot(m_N.T, Phi_x)#予測分布の平均m_N^T\phi(x)
            s2_N=1.0/beta+dot(Phi_x.T,solve(S_N_inv,Phi_x))#予測分布の分散1/\beta+\phi(x)^TS_N\Phi(x)
            return (mu,s2_N)
        
        return m_N,S_N_inv,func#(3.53),(3.54),(3.58)のxに対する平均、分散を返す関数
    
    xmin=-0.05
    xmax=1.05
    ymin=-1.5
    ymax=1.5
    
    
    clf()
    axis([xmin,xmax,ymin,ymax])
    title("Fig 3.8 (%d sample%s)" % (N,'s' if N > 1 else ''))
    
    x_=arange(xmin,xmax,0.01)
    plot(x_,sin(x_*pi*2),color='gray')#元のデータ

    m_N,S_N_inv,f=predictive_dist_func(ALPHA,BETA)#ハイパーパラメータに対応する平均、分散、予測式
    
    y_h=[]
    y_m=[]
    y_l=[]
    #yの予測の平均(y_m)、平均±標準偏差(y_h,y_l)作成
    for mu,s2 in map(f,x_):#mu,s2:各x_の要素に対応する予測の平均と分散
        s=sqrt(s2)
        y_m.append(mu)
        y_h.append(mu+s)
        y_l.append(mu-s)
    #図3.8描画
    fill_between(x_,y_h,y_l,color='#cccccc')
    plot(x_,y_m,color='#000000')
    scatter(xs,ts,color='r', marker='o')#入力と目標値のペア
    show()
    clf()
    
    x_=arange(xmin,xmax,0.01)
    #図3.9描画
    plot(x_,sin(x_*pi*2),color='gray')
    #推定したwでグラフを描画
    for i in range(5):
        w=multivariate_normal(m_N,inv(S_N_inv),1).T#wを作成(平均m_N、分散S_Nガウス分布より)
        y=lambda x:dot(w.T,Phi(x))[0]#推定値を算出する関数
        plot(x_,y(x_),color='#cccccc')
    scatter(xs,ts,color='r',marker='o')
    show()
 
    
def main():
    #xと対応するtを作成(t_n=sin(x_n*2*\pi)+ガウスノイズ)
    xs=arange(0,1.01,0.02)#x_n
    ts=sin(xs*pi*2)+normal(loc=0.0,scale=0.1,size=size(xs))#t_n
    
    #n:範囲、k:個数
    #0~n-1からk個選んでソートしたものを返す。抜き出すデータ点
    def randidx(n,k):
        r = list(range(n))
        shuffle(r)
        return sort(r[0:k])
    #それぞれの訓練データ数(1,2,5,20)で処理
    for k in(1,2,5,20):
        indices=randidx(size(xs),k)
        sub(xs[indices],ts[indices])#使用する(x_n,t_n)
        
if __name__=='__main__':
    main()