prml4.tex

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 第4章
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\setcounter{chapter}{3}
\chapter{「線形識別モデル」のための数学}
\section{クラス分類問題}
クラス分類とは, 与えられた入力空間を相異なる$K$個の空間に分割し, それぞれの空間をクラス$C_k$とラベルをつけることである.
訓練データ$x$が与えられたときに, 推論段階と決定段階を経て各クラスに割り当てる.
$$
\text{訓練データ}x
\rightarrow \text{モデル$p(C_k\,|\,x)$を作る}
\rightarrow \text{事後確率を使ってクラスに割り当てる}
$$

訓練データからどの情報を使って分類するかによって三つの方法がある：
\begin{itemize}
\item 生成モデル（generative model）：
$p(x\,|\,C_k)$を$C_k$ごとに決め, $p(C_k)$も決める. そうすると同時分布$p(x, C_k)$が分かり,
$$
p(C_k\,|\,x)=\frac{p(x\,|\,C_k)\,p(C_k)}{p(x)}
$$
で事後確率を求める.

$p(x)$は$p(x)=\sum p(x\,|\,C_k)\,p(C_k)$で求められる.
$p(x\,|\,C_k)$があると自分でさいころを振って各$C_k$に対して$x$を作ることができるという点で,
生成モデルという.

\item 識別モデル（discriminative model）：
$p(x\,|\,C_k)$を求めずにいきなり事後確率$p(C_k\,|\,x)$を決める推論問題を解く.
決定理論を使って$x$をあるクラスに割り当てる.

\item 識別関数（discriminant function）：
確率モデルを考えずに入力関数によって定まる識別関数$f\colon x \mapsto k$を作る.
\end{itemize}

\section{行列の微分の復習}
$A=(a_{ij})$と書いた.
$$
(AB)_{ij}=\sum_k a_{ik} b_{kj}, \quad \tr(A)=\sum_i a_{ii}, \quad \trans{A}=(a_{ji})
$$
などを思い出しておく.
さて$A$, $B$を適当な行列として
$$
\dif{A}\tr(AB)=\trans{B}
$$
なぜなら,
$$
\left(\dif{A}\tr(AB)\right)_{ij}=\dif{a_{ij}}\sum_{s,t} a_{st} b_{ts} = b_{ji}.
$$
ここで
$\partial a_{st}/ \partial a_{ij}= \delta_{is} \delta_{jt}$を使った. つまり添え字$s$, $t$が走るときに, $s=i$, $t=j$のときのみが生き残るというわけである.
慣れるためにもう一つやっておこう.
$$
\dif{A}\tr(AB\trans{A})=A(B+\trans{B}).
$$
なぜなら,
\begin{eqnarray*}
\dif{a_{ij}} \tr(AB\trans{A})
 &=& \dif{a_{ij}} \sum_{s,t,u} a_{st} b_{tu} a_{su}
 = \sum_{s,t,u} b_{tu}  \dif{a_{ij}} (a_{st} a_{su})\\
 &=& \sum_{s,t,u} b_{tu} \left(\delta_{is} \delta_{jt} a_{su} + a_{st} \delta_{is} \delta_{ju}\right)
 = \sum_{u} b_{ju} a_{iu} + \sum_{t} b_{tj} a_{it}\\
 &=& \sum_{u} a_{iu} b_{ju} + \sum_{t} a_{it} b_{tj}
 = (A\trans{B})_{ij} + (AB)_{ij}
 = \left(A(B + \trans{B})\right)_{ij}.
\end{eqnarray*}
\vspace{0pt}

\section{多クラス}
$K$個の線形関数を使った$K$クラス識別を考える.
$$
y_k(x)=\trans{w_k}x+w_{k0}.
$$
ここで$w_k$は重みベクトル, $w_{k0}$はバイアスパラメータでスカラー, $x$が分類したい入力パラメータでベクトルである.
クラス分類を次の方法で定義する:
$x$に対して, ある$k$が存在し, 全ての$j\neq k$に対して$y_k(x) > y_j(x)$であるとき$x$はクラス$C_k$に割り当てるとする.
これはwell-definedである.
つまり
\begin{itemize}
\item  （一意性）$x$が二つの異なるクラス$C_k$に$C_k'$に属することはない. なぜならそういう$k$, $k'$があったとすると$y_k(x) > y_k'(x) > y_k(x)$となり矛盾するから.
\item （存在性）$x$が与えられたとき$\{y_k(x)\}$の最大値$m$を与える$k_0$がその候補である. もしも
$m=y_k(x)$となる$k$が複数個存在（$k_1$, $k_2$）したとすると,
クラス分類はできないが, そういう$x$の集合は$\{x\,|\,y_{k_1}(x)=y_{k_2}(x)\}$の
部分集合となり, 通常次元が落ちる. つまり無視できるぐらいしかない.
\end{itemize}
上記で分類されたクラス$C_k$に属する空間は凸領域となる.
すなわち
$x$, $x'$を$C_K$の点とすると, 任意の$\lambda \in [0, 1]$に対して$x''=\lambda x + (1-\lambda) x'$も$C_k$に属する.

なぜなら$x$, $x' \in C_k$より任意の$j\neq k$に対して$y_k(x) > y_j(x)$, $y_k(x') > y_j(x')$.
$y_k(x)$は$x$について線形なので$\lambda \geq 0$, $1-\lambda \geq 0$より
$$
y_k(x'') = \lambda y_k(x) + (1-\lambda) y_k(x') > \lambda y_j(x) + (1-\lambda) y_j(x') = y_j(x'')
$$
が成り立つからである.

領域内の任意のループを連続的に1点につぶすことが出来る様な領域を単連結（simply connected）という.
凸領域は単連結である. つまりその領域の中に空洞は無い. 任意の凸領域の2点を結ぶ線分が凸領域に入ることから直感的には明らかであろう.

単連結であることを簡単に示しておこう：$X$を凸領域, $S^1=\{(x,y)\,|\,x^2+y^2=1\}$を単位円とする.
$$
f:S^1 \rightarrow X
$$
を$S^1$からXへの連続関数とする. 任意のループは$f(S^1)$で表される.
$t \in [0, 1]$, $\lambda \in [0, 1]$に対して
$$
f_\lambda(t) := \lambda f(t) + (1-\lambda)f(0)
$$
とすると$X$の凸性から$f_\lambda(t) \in X$.
つまり$f_\lambda(S^1)$は$X$内のループ.
$f_0(t)=f(0)$は$X$のある1点. $f_1(t)=f(t)$は元のループだから, これはループ$f(S^1)$を$X$の中で連続的に一点$f(0)$につぶすことが出来ることを示している. つまり$X$は単連結.

\section{分類における最小二乗}
前節では重みベクトル$w_{k0}$を別扱いしたが, $\tilde{w}_k:=\trans{(w_{k0}, \trans{w_k})}$,
$\tilde{x}:=\trans{(1,\trans{x})}$と1次元増やすと$y_k(x)=\trans{\tilde{w}}\tilde{x}$と書ける.
面倒なので$\tilde{x}$を$x$と置き換えてしまおう.

さらにまとめて
$y(x)=\trans{W}x$としよう. $x$, $y$はベクトル, $W$は行列である.

二乗誤差関数
$$
E_D(W):=\half\tr\left(\trans{(XW-T)}(XW-T)\right)
$$
を最小化する$W$を求めよう.
\begin{eqnarray*}
\dif{w_{ij}}E_D(W)
 &=& \half\dif{w_{ij}}\sum_{s,t} \left((XW-T)_{st}\right)^2
 = \sum_{s,t} (XW-T)_{st} \dif{w_{ij}}(XW-T)_{st}\\
 &=& \sum_{s,t} (XW-T)_{st} \dif{w_{ij}}\left(\sum_u x_{su}w_{ut}\right)
 = \sum_{s,t,u} (XW-T)_{st} x_{su}\delta_{iu}\delta_{jt}\\
 &=& \sum_s (XW-T)_{sj} x_{si}
 = \sum_s (\trans{X})_{is} (XW-T)_{sj}
 = \left(\trans{X}(XW-T)\right)_{ij}.
\end{eqnarray*}
よって
$$
\dif{W}E_D(W)=\trans{X}(XW-T).
$$
$=0$とおいて
$\trans{X}XW = \trans{X}T$より
$$
W=(\trans{X}X)^{-1}\trans{X}T.
$$
\vspace{0pt}

\section{フィッシャーの線形判別}
まず$D$次元のベクトル$x$の入力に対して$y=\trans{w}x$で1次元に射影する.
$y \ge w_0$なら$C_1$, そうでないなら$C_2$に分類する.
$C_1$の点が$N_1$個, $C_2$の点が$N_2$個とする.
$C_i$の点の平均は
$$
\bmm_i:=\frac{1}{N_i}\sum_{n \in C_i} x_n.
$$
$m_i=\trans{w}\bm{m_i}$として, $|w|^2=\sum_i w_i^2=1$の制約下で
$$
m_2 - m_1=\trans{w}(\bmm_2-\bmm_1)
$$
を最大化してみよう.
ラグランジュの未定乗数法を用いて
$$
f(w):=\trans{w}(\bmm_2-\bmm_1)+\lambda(1-|w|^2)
$$
とおくと
$$
\diff{w}{f}=\bmm_2-\bmm_1-2\lambda w = 0.
$$
よって
$$
w=\frac{1}{2\lambda}(\bmm_2-\bmm_1) \propto (\bmm_2-\bmm_1).
$$
$$
\diff{\lambda}{f}=1-|w|^2=0
$$
より$|w|=1$.
ただしこの手法ではそれぞれのクラスの重心$\bmm_1$と$\bmm_2$とだけで$w$の向きが決まってしまい,
場合によっては二つのクラスの射影が大きく重なってうまく分離できないことがある.
そこでクラス間の重なりを最小にするように分散も加味してみる.

クラス$C_k$から射影されたデータのクラス内の分散を
$$
y_n:=\trans{w}x_n, \quad s_k^2:=\sum_{n\in C_k}(y_n-m_k)^2
$$
で定義し, 全データに対する分散を$s_1^2+s_2^2$とする.
フィッシャーの判別基準は
$$
J(w):=\frac{(m_2-m_1)^2}{s_1^2+s_2^2}
$$
で定義される. この定義を書き直してみよう.
\begin{eqnarray*}&&
S_B:=(\bmm_2-\bmm_1)\trans{(\bmm_2-\bmm_1)},
\\&&
S_W:=\sum_{n \in C_1}(x_n-\bmm_1)\trans{(x_n-\bmm_1)}+\sum_{n \in C_2}(x_n-\bmm_2)\trans{(x_n-\bmm_2)}
\end{eqnarray*}
とする.
$S_B$をクラス間共分散行列, $S_W$を総クラス内共分散行列という.
\begin{eqnarray*}&&
\trans{w}S_B w
 = \trans{w}(\bmm_2-\bmm_1)\trans{(\bmm_2-\bmm_1)}w
 = (m_2-m_1)^2,
\\&&
\trans{w}S_W w
 = \trans{w}\sum_{n \in C_1}(x_n-\bmm_1)\trans{(x_n-\bmm_1)}w+\trans{w}\sum_{n \in C_2}(x_n-\bmm_2)\trans{(x_n-\bmm_2)}w
\\&&\phantom{\trans{w}S_W w }
 = \sum_{n\in C_1} (y_n-m_1)^2+\sum_{n\in C_2}(y_n-m_2)^2
\end{eqnarray*}
より
$$
J(w)=\frac{\trans{w}S_B w}{\trans{w} S_W w}.
$$
これが最大となる$w$の値を求めてみよう.
大きさはどうでもよくて向きが重要である.
$$
\dif{w}J(w)=\left(2 (S_B w)(\trans{w} S_W w) - 2 (\trans{w} S_B w) (S_W w)\right)/(\trans{w}S_W w)^2=0.
$$
よって
$$
(\trans{w}S_B w)S_w w = (\trans{w} S_W w)S_B w.
$$
$S_B w = (\bmm_2-\bmm_1)\left(\trans{(\bmm_2-\bmm_1)}w\right)\propto (\bmm_2-\bmm_1)$だから
$$
w \propto S_W^{-1}S_B w \propto S_W^{-1}(\bmm_2-\bmm_1)
$$
のときに$J(w)$が最大となる. これをフィッシャーの線形判別（linear discriminant）という.

\section{最小二乗との関連}
\begin{itemize}
\item 最小二乗法：目的変数の値の集合にできるだけ近いように
\item フィッシャーの判別基準：クラスの分離を最大化するように
\end{itemize}
2クラスの分類のときは最小二乗の特別な場合がフィッシャーの判別基準であることをみる.
フィッシャーの判別基準が, 最小二乗と関係があることが分かるとそちらの議論が使えていろいろ便利なことがある.

クラス$C_i$に属するパターンの個数を$N_i$として全体を$N:=N_1+N_2$とする.
クラス$C_1$に対する目的変数値を$N/N_1$, クラス$C_2$に対する目的変数値を$-N/N_2$とする.

この条件下で二乗和誤差
$$E:=\half\sum_{n=1}^N(\trans{w}x_n+w_0-t_n)^2$$
を最大化してみよう.
$$
\diff{w_0}{E}=\sum(\trans{w}x_n+w_0-t_n)=0
$$
より$\bmm:=(1/N)\sum x_n$とおくと$N \trans{w}m+N w_0 - \sum t_n=0$.
$$\sum t_n=N_1(N/N_1) + N_2(-N/N_2)=0$$
より$w_0=-\trans{w}\bmm$. また
$$
\sum(\trans{w}x_n)x_n=\sum(\trans{x_n}w)x_n=\sum(x_n \trans{x_n})w.
$$
$$
\sum w_0 x_n = N w_0 \bmm = -N(\trans{w}\bmm)\bmm=-N(\bmm \trans{\bmm})w.
$$
$$
\sum t_n x_n = \sum_{n\in C_1} t_n x_n + \sum_{n\in C_2} t_n x_n
= N/N_1(N_1 \bmm_1)+(-N/N_2)(N_2 \bmm_2)=N(\bmm_1-\bmm_2).
$$
よって
$$
\dif{w}E=\sum\left(\trans{w}x_n + w_0-t_n\right)x_n=0
$$
を使うと
$$
\sum(x_n \trans{x_n})w=N(\bmm \trans{\bmm})w+N(\bmm_1-\bmm_2).
$$
これらの式を使って$S_w$を計算する.
\begin{eqnarray*}
S_W
 &=& \sum_{n\in C_1} x_n\trans{x_n}-2\sum_{C_1} x_n \trans{\bmm_1}+\sum_{C_1} \bmm_1\trans{\bmm_1} + \sum_{C_2} x_n\trans{x_n} - 2\sum_{C_2} x_n\trans{\bmm_2}\\
 &+& \sum_{C_2} \bmm_2 \trans{\bmm_2}
 = \sum x_n\trans{x_n}-N_1 \bmm_1 \trans{\bmm_1}-N_2 \bmm_2\trans{\bmm_2}\\
 &=&N(\bmm \trans{\bmm})w+N(\bmm_1-\bmm_2)-N_1 \bmm_1 \trans{\bmm_1}-N_2 \bmm_2\trans{\bmm_2}.
\end{eqnarray*}
よって
\begin{eqnarray*}
\left(S_W+\frac{N_1 N_2}{N}S_B\right)w&=&N(\bmm_1-\bmm_2)\\
                           &+&\Bigl(N\bmm \trans{\bmm}-N_1\bmm_1 \trans{\bmm_1}-N_2 \bmm_2 \trans{\bmm_2}\\
                     &+&\frac{N_1 N_2}{N} (\bmm_1 -\bmm_2)\trans{(\bmm_1-\bmm_2)}\Bigr)w.
\end{eqnarray*}

$()$内が$0$であることを示す．
\begin{eqnarray*}&&
()
 = \frac{1}{N}\outp{\left(N_1 \bmm_1+N_2 \bmm_2\right)}-N_1 \outp{\bmm_1}-N_2 \outp{\bmm_2}
\\&&\phantom{()={}}
   + \frac{N_1 N_2}{N}\left(\bmm_1 \trans{\bmm_2}+\bmm_2 \trans{\bmm_2}\right)
\\&&\phantom{()}
 = \left(\frac{N_1^2}{N}-N_1+\frac{N_1 N_2}{N}\right)\outp{\bmm_1}+\left(\frac{2}{N} N_1 N_2 - \frac{2}{N}N_1 N_2\right)\bmm_1\trans{\bmm_2}
\\&&\phantom{()={}}
   + \left(\frac{N_2^2}{N}-N_2+\frac{N_1 N_2}{N}\right)\outp{\bmm_2},
\\&&
\frac{N_1^2}{N}-N_1+\frac{N_1 N_2}{N}=\frac{N_1}{N}(N_1-N+N_2)=0,
\\&&
\frac{N_2^2}{N}-N_2+\frac{N_1 N_2}{N}=\frac{N_2}{N}(N_2-N+N_1)=0.
\end{eqnarray*}
よって
$$
\left(S_W+\frac{N_1 N_2}{N}S_B\right)w=N(\bmm_1-\bmm_2).
$$
$S_Bw\propto(\bmm_2 -\bmm_1)$なので$w \propto S_W^{-1}(\bmm_2 - \bmm_1)$.

\section{確率的生成モデル}
分類を確率的な視点から見る.
生成的アプローチ
\begin{itemize}
\item $p(x\,|\,C_k)$：モデル化されたクラスの条件付き確率密度
\item $p(C_k)$：クラスの事前確率
\end{itemize}
$$
p(C_1\,|\,x)=\frac{p(x\,|\,C_1)\,p(C_1)}{p(x\,|\,C_1)\,p(C_1)+P(x\,|\,C_2)\,p(C_2)}
$$
とする. ロジスティックシグモイド関数を
$$
\sigma(a):=\frac{1}{1+\exp(-a)}
$$
と定義し,
$$
a:=\log \frac{p(x\,|\,C_1)\,p(C_1)}{p(x\,|\,C_2)\,p(C_2)}とすると
$$
$$
\sigma(a)=\frac{1}{1+\frac{p(x\,|\,C_2)\,p(C_2)}{p(x\,|\,C_1)\,p(C_1)}}=p(C_1\,|\,x).
$$

ロジスティックシグモイド関数の関数の性質：
$$
\sigma(-a)=\frac{1}{1+e^a}=1-\frac{e^a}{1+e^a}=1-\frac{1}{1+e^{-a}}=1-\sigma(a).
$$
$$
\sigma(a)=\frac{e^a}{e^a+1}
$$
より$e^a(\sigma(a)-1)=-\sigma(a)$. よって
$$
a=\log \frac{\sigma(a)}{1-\sigma(a)}.
$$
この関数をロジット関数という.

$K>2$クラスの場合, $a_k=\log(p(x\,|\,C_k)\,p(C_k))$より
$$
p(C_k\,|\,x)=\frac{p(x\,|\,C_k)\,p(C_k)}{\sum_j p(x\,|\,C_j)\,p(C_j)}=\frac{\exp(a_k)}{\sum_j \exp(a_j)}
$$
となる. この関数は正規化指数関数, あるいはソフトマックス関数という.

\section{連続値入力}
条件付き確率密度がガウス分布, そのガウス分布の共分散行列（$\Sigma=A$）がすべてのクラスで共通と仮定する.

$$
p(x\,|\,C_k)=\frac{1}{(2\pi)^{D/2}}\frac{1}{|A|^{1/2}}\exp\left(-\half\trans{(x-\mu_k)}A^{-1}(x-\mu_k)\right).
$$
\begin{eqnarray*}
a&=& \log \frac{p(x\,|\,C_1)\,p(C_1)}{p(x\,|\,C_2)\,p(C_2)}\\
 &=& \log \frac{p(C_1)}{p(C_2)}-\half\trans{(x-\mu_1)}A^{-1}(x-\mu_1)+\half\trans{(x-\mu_2)}A^{-1}(x-\mu_2)\\
 &=& \log \frac{p(C_1)}{p(C_2)}-\half\quads{A^{-1}}{\mu_1}+\half\quads{A^{-1}}{\mu_2}+\trans{(\mu_1-\mu_2)}Ax.
\end{eqnarray*}
よって
$$
w_0:=-\half\quads{A^{-1}}{\mu_1}+\half\quads{A^{-1}}{\mu_2}+\log \frac{p(C_1)}{p(C_2)},
$$
$$
w:=A^{-1}(\mu_1-\mu_2)
$$
とおくと$p(C_1\,|\,x)=\sigma(\trans{w}x+w_0)$. つまりロジスティックシグモイド関数の中は$x$について線形.

$K$クラスの場合, 上で定義した$a_k$を用いると
\begin{eqnarray*}
a_k
 &=&\log p(C_k)-\half\quads{A^{-1}}{\mu_k}+\trans{\mu_k}A^{-1}x-\half\trans{x}A^{-1}x+\text{const.}\\
 &=& a_k'-\half\trans{x}A^{-1}x+\text{const.}
\end{eqnarray*}
ここで$a_k':=\trans{w_k}x+w_{k0}$, $w_k:=A^{-1}\mu_k$, $w_{k0}:=-(1/2)\quads{A^{-1}}{\mu_k}+\log p(C_k)$.

（注）PRML (4.63)の$a_k$の定義だと$x$の2次の項が残るため, 式(4.68)を出すにはそれを除かなければならない.

よって
\begin{eqnarray*}&&
p(C_k\,|\,x)
=\frac{\exp(a_k)}{\sum_j \exp(a_j)}
=\frac{\exp(a_k')\exp\left(-(1/2)\quads{A^{-1}}{x}+\text{const.}\right)}{\sum_j \exp(a_j')\exp\left(-(1/2)\quads{A^{-1}}{x}+\text{const.}\right)}
=\frac{\exp(a_k')}{\sum_j \exp(a_j')}.
\end{eqnarray*}
\vspace{0pt}

\section{最尤解}
条件付き確率分布がガウス分布, それらが共通の共分散行列を持つと仮定する.

2クラスの場合を考える. 
データ集合$\{(x_n, t_n)\}$, $n=1, \ldots, N$. $t_n=1$はクラス$C_1$,
$t_n=0$はクラス$C_2$とする\footnote{PRMLでは「データ集合$\{x_n, t_n\}$」と
書かれているが, $x_n$と$t_n$は対であるため括弧で囲った.}.
さらに$p(C_1)=p$, $p(C_2)=1-p$という事前確率を割り当てる. $N_i$をクラス$C_i$のデータの個数, $N=N_1+N_2$を総数とする.
\begin{eqnarray*}
&& p(x_n, C_1)=p(C_1)\,p(x_n\,|\,C_1)=p \,\calN(x_n\,|\,\mu_1, A),
\\
&& p(x_n, C_2)=p(C_2)\,p(x_n\,|\,C_2)=p (1-p) \,\calN(x_n\,|\,\mu_2, A).
\end{eqnarray*}
尤度関数は
$$
p(t,X\,|\,p, \mu_1, \mu_2, A):= \prod^N\left(p \,\calN(x_n\,|\,\mu_1, A)\right)^{t_n} \left((1-p)\,\calN(x_n\,|\,\mu_2,A)\right)^{1-t_n}.
$$
このうち$p$に関する部分の対数は
$$
\sum \left(t_n \log p + (1-t_n) \log (1-p)\right).
$$
$p$で微分して$0$とおく.
$$
\frac{1}{p}\sum t_n - \frac{1}{1-p}\sum (1-t_n)=\frac{1}{p} N_1 - \frac{1}{1-p} N_2=\frac{(1-p)N_1-pN_2}{p(1-p)}=0.
$$
よって$p=N_1/(N_1+N_2)=N_1/N$.
つまり$p$に関する最尤推定は$C_1$内の個数になる.

$K$クラスのときを考えてみよう.
$\sum p_i=1$. 尤度関数は
$$
p(t,X\,|\,p_1, \ldots, p_K, \mu_1, \ldots, \mu_K, A)=\prod_{n=1}^N \prod_{i=1}^K p_i \,\calN(x_n\,|\,\mu_i,A))^{t_{ni}}.
$$
この対数に未定乗数法の$\lambda (\sum p_i - 1)$の項を加え, $p_i$に関する部分を抜き出すと
$$
\sum_n t_{ni} \log p_i+\lambda p_i.
$$
$p_i$で微分して$0$とおくと$\sum_n (t_{ni}/p_i) + \lambda=0$.
よって
$$
-p_i \lambda = \sum_n t_{ni}=N_i
$$
また
$-\sum_i p_i \lambda = -\lambda = \sum_i N_i = N$より$p_i=-N_i/\lambda=N_i/N$.

さて2クラスの問題に戻って$\mu_i$について最大化してみよう. $\mu_1$についての部分は
$$
\sum t_n \log \calN(x_n\,|\,\mu_1,A)=-\half\sum t_n \quads{A^{-1}}{(x_n-\mu_1)} + \text{const.}
$$
$\mu_1$で微分して$0$とおくと
$$
\sum t_n A^{-1}(x_n-\mu_1)=A^{-1} \left(\sum t_n x_n - \mu_1 \sum t_n\right)=A^{-1} \left(\sum t_n x_n - \mu_1 N_1\right)=0.
$$
よって
$$\mu_1=\frac{1}{N_1}\sum t_n x_n.
$$
$\mu_2$については$\sum (1-t_n) \log \calN(x_n\,|\,\mu_2,A)$を考えて
$$
\mu_2=\frac{1}{N_2}\sum (1-t_n)x_n.
$$
最後に$A$に関する最尤解を求める. $A$に関する部分の対数は
\begin{eqnarray*}
 && -\half\sum_{n=1}^N\Bigl(t_n \log |A|+t_n \trans{(x_n-\mu_1)}A^{-1}(x_n-\mu_1)\\
 &&+(1-t_n) \log |A|+(1-t_n)\trans{(x_n-\mu_2)}A^{-1}(x_n-\mu_2)\Bigr)\\
 &=& -\frac{N}{2}\log |A|\\
 &&{}-\half \tr\left(\sum_{n=1}^N\left(t_n A^{-1}(x_n -\mu_1)\trans{(x_n-\mu_1)}+(1-t_n)A^{-1}(x_n-\mu_2)\trans{(x_n-\mu_2)}\right)\right)\\
 &=& -\frac{N}{2}\log |A|\\
 &&{}-\half \tr\left(A^{-1}\left(\sum_{n\in C_1}(x_n-\mu_1)\trans{(x_n-\mu_1)}+\sum_{n\in C_2}(x_n-\mu_2)\trans{(x_n-\mu_2)}\right)\right)\\
 &=& -\frac{N}{2}\log |A| -\half \tr\left(A^{-1}(N_1 S_1 + N_2 S_2)\right)
 = -\frac{N}{2}\log |A| -\frac{N}{2} \tr(A^{-1}S)
\end{eqnarray*}
ここで最後の式変形に
$$
S_i:=\frac{1}{N_i}\sum_{n\in C_i}(x_n-\mu_i)\trans{(x_n-\mu_i)},
\hseq
S:=\frac{N_1}{N}S_1+\frac{N_2}{N} S_2
$$
を用いた. これを$A$で微分する.
2章で示した行列式の対数の微分の公式(2)：式(\ref{diff_log_mat})
$$
\dif{A}\log|A|=\trans{(A^{-1})}
$$
と式(\ref{diff_tr_invA_B})：
$$
\dif{A}\tr(A^{-1}B)=-\trans{(A^{-1}BA^{-1})}
$$
を使うと
$$
-\frac{N}{2}\left(\trans{(A^{-1})}-\trans{(A^{-1}SA^{-1})}\right)=0.
$$
よって$A=S$となる. これは2クラスの各クラスの共分散行列の重みつき平均である.
またフィッシャーの判別基準で求めた総クラス内共分散行列$S_W$を$N$で割ったものに等しいことにも注意する.

\section{ロジスティック回帰}
2クラス分類問題において, ある程度一般的な仮定の下で$C_1$の事後確率は
$$
p(C_1\,|\,\phi)=y(\phi)=\sigma(\trans{w}\phi)
$$
とかけた. もちろん$p(C_2\,|\,\phi)=1-p(C_1\,|\,\phi)$である.
この式の導出に使った仮定を忘れ, これを出発点としこの形の関数を使うモデルをロジスティック回帰（logistic regression）という.
このモデルにおけるパラメータを最尤法で求める.
$$
\sigma(x)=\frac{1}{1+e^{-x}}
$$
としたとき
\pagebreak
$$
\sigma'(x)
=-\frac{-e^{-x}}{(1+e^{-x})^2}
=\frac{e^{-x}}{(1+e^{-x})^2}
=\frac{1}{1+e^{-x}}\left(1-\frac{1}{1+e^{-x}}\right)
=\sigma(x)(1-\sigma(x)).
$$
データ集合$\{(\phi_n, t_n)\}$,
$t_n\in \{0,1\}$, $\phi_n=\phi(x_n)$, $n=1, \ldots, N$,
$t=\trans{(t_1, \ldots, t_N)}$, $y_n=p(C_1\,|\,\phi_n)=\sigma(a_n)$, $a_n=\trans{w}\phi_n$とする.
尤度関数は
$$
p(t\,|\,w)=\prod_{n=1}^N y_n^{t_n}(1-y_n)^{1-t_n}.
$$
誤差関数は
$$
E(w)=-\log p(t\,|\,w)=-\sum\left(t_n \log y_n + (1-t_n) \log (1-y_n)\right).
$$
$w$で微分してみよう.
まず
$$
\dif{w}y_n=\sigma'(a_n) \phi_n=y_n(1-y_n) \phi_n.
$$
よって
\begin{eqnarray*}
\dif{w}E &=&-\sum\left(t_n (1-y_n)\phi_n + (1-t_n)(-y_n)\phi_n\right)\\
 &=& \sum(-t_n+t_n y_n+y_n-y_n t_n)\phi
 = \sum(y_n-t_n)\phi_n.
\end{eqnarray*}
$y_n-t_n$は目的値とモデルの予測値との誤差なので,
PRML~3.1.1で扱われた線形回帰モデルのときと同じ形になる.

\section{反復再重み付け最小二乗}
いわゆるニュートン・ラフラソン法は関数$f(x)$の零点を求める方法である：
零点の近似解$x_n$が与えられたときにより近い値$x_{n+1}$を見つける.
$x_n$における接線の方程式$
f'(x_n)(x-x_n)+f(x_n)=0
$
の解を$x_{n+1}$とすると
$
x_{n+1}=x_n - (f'(x_n))^{-1} f(x_n)
$であるが，ここで$n\to\infty$として
得られる$x=\lim_{n\to\infty} x_n$が$f(x)$の零点である.

さて，関数$E(w)$を最小化するためのベクトル$w$を与える更新式を
考えてみると, 最小化を与える$w$は$\dif{w}E(w)$の零点.
$\dif{w}E(w)$の零点を求める問題にニュートン・ラフラソン法を適用する.
$w$を古い値, $w'$を新しい値, $H(w)$を$E(w)$のヘッシアンとする.
$f \longleftrightarrow \dif{w}E(w)$, $f' \longleftrightarrow H(w)$という対応により
$$
w'=w-H(w)^{-1}\dif{w}E(w).
$$
この式を線形回帰モデルに適用してみる.
$$
E(w)=\half\sum\left(t_n-\trans{w}\phi(x_n)\right)^2, \phi_n=\phi(x_n)
$$
を$w$で微分して$\Phi=\trans{(\phi_1, \ldots)}$とおくと
$$
\dif{w}E(w)=\sum(t_n-\trans{w}\phi_n)\phi_n=\sum \phi_n \trans{\phi_n}w-\sum \phi_n t_n=\PhiT \Phi w-\PhiT t.
$$
$$
H=H(w)=\ddif{w_i}{w_j}E(w)=\PhiT \Phi.
$$
更新式に代入すると
$$
w'=w-(\PhiT \Phi)^{-1}(\PhiT \Phi w-\PhiT t)=(\PhiT \Phi)^{-1}\PhiT t.
$$
これは最小二乗解である. つまり1回の更新で厳密解に到達した. 線形回帰モデルは$w$に関して2次なので$\dif{w}E(w)$に関しては1次だからである.
次にロジスティック回帰に適用してみる.
$$
E(w)=-\sum\left(t_n \log y_n + (1-t_n)\log (1-y_n)\right), y_n=\sigma(a_n)=\sigma(\trans{w}\phi_n).
$$
$$
\dif{w}E(w)=\sum (y_n-t_n)\phi_n.
$$
$y_n'=y_n(1-y_n)\phi_n$だったので
$$
H=\sum \phi_n y_n(1-y_n)\trans{\phi_n}=\PhiT R\Phi.
$$
ここで$R=\diag(R_n)=\diag(y_n(1-y_n))$.

$H>0$を確認する.
任意の縦ベクトル$u$に対して$v=\Phi u$とおくと$v\ne 0$ならば
$$
\trans{u}H u=\trans{v}R v=\sum y_n(1-y_n)v_n^2 > 0.
$$
最後の不等号では$0<y_n<1$を用いた. ヘッシアンが正定値であることが分かったので交差エントロピー誤差関数は唯一の最小解を持つ.

$w$の更新式を見てみよう.
\begin{eqnarray*}
w'
 &=& w-(\PhiT R\Phi)^{-1}\PhiT(y-t)
 = (\PhiT R\Phi)^{-1}(\PhiT R\Phi w - \PhiT(y-t))\\
 &=& (\PhiT R\Phi)^{-1}\PhiT R(\Phi w - R^{-1}(y-t))
 = (\PhiT R\Phi)^{-1}\PhiT Rz.
\end{eqnarray*}
ここで$z=\Phi w - R^{-1}(y-t)$である.
$R$は$y_n$つまり$w$に依存しているので正規方程式は更新式ごとに計算し直す必要がある.
反復再重み付き最小二乗法（IRLS: iterative reweighted least squares method）という.

$t=1$をクラス$C_1$, $t=0$をクラス$C_2$に割り当てて, それぞれの確率は$y$, $1-y$だから
$$
\EE[t]=y=\sigma(x).
$$
$t^2=t$だから
$$
\var[t]=\EE[t^2]-\EE[t]^2=\EE[t]-\EE[t]^2=y-y^2=y(1-y).
$$
つまり重み付け対角行列$R$の対角成分は分散である.

IRLSを線形近似の解として解釈することも出来る. すなわち$a=\trans{w}\phi$, $y=\sigma(a)$という関係を通じて
$a$を$y$の関数とみなし,
$a_n$を目標値$t_n$の変数とみなして近次解$y_n=\sigma(\trans{w_{\text {old}}}\phi)$のまわりで一次近似を行うと
\begin{eqnarray*}
a_n &\approx& a_n(y_n)+\left.\dif{y_n} a_n \right|_{t_n=y_n}(t_n-y_n)
 = \trans{w_{\text {old}}}\phi + \frac{1}{y_n(1-y_n)}(t_n-y_n)\\
 &=& \trans{w_{\text {old}}}\phi - \frac{y_n-t_n}{y_n(1-y_n)}
 = z_n.
\end{eqnarray*}
つまり$z_n$は線形近似したときの目標変数値と解釈できる.

\section{Jensenの不等式}
実数上の実数値関数$f(x)$が凸関数であるとする.
すなわち任意の$x$, $y$, $0 \le t \le 1$に対して
$$
tf(x) + (1-t)f(y) \ge f(tx+(1-t)y)
$$
である.

$p_1, \ldots, p_n$を足して1になる非負の数, すなわち$\sum_{i=1}^n p_i=1$, $p_i \ge 0$とする.

このとき$n$個の任意の実数$x_1, \ldots, x_n$に対して
$$
\sum_{i=1}^n p_i f(x_i) \ge f\left(\sum_{i=1}^n p_i x_i\right).
$$
これをJensenの不等式という.

証明は数学的帰納法を使う.
$n=1$のときは自明. $n$のとき成り立つとし,
$$
\sum_{i=1}^{n+1} p_i=1
$$
とする.
$q=\sum_{i=1}^n p_i$とおくと$q + p_{n+1}=1$.
$q=0$のときは$p_{n+1}=1$となり上記不等式は自明になりたつ.
\pagebreak
$q \ne 0$のときは
$
\sum_{i=1}^n (p_i/q) = 1
$に注意して計算すると,
\begin{eqnarray*}
\sum_{i=1}^{n+1} p_i f(x_i) &=& q \sum_{i=1}^n (p_i/q) f(x_i) + p_{n+1} f(x_{n+1})\\
 &&\text{（帰納法の仮定を用いて）}\\
 &\ge& q f\left(\sum_{i=1}^n (p_i/q) x_i\right) + p_{n+1} f(x_{n+1})\\
 &&\text{（$f$が凸関数であることを用いて）}\\
 &\ge& f\left(q \left(\sum_{i=1}^n (p_i/q) x_i\right) + p_{n+1}x_{n+1}\right)
 = f\left(\sum_{i=1}^{n+1} p_i x_i\right).
\end{eqnarray*}
\vspace{0pt}

\section{多クラスロジスティック回帰}\label{takurasu}
多クラス分類の事後確率を
$$
p(C_k\,|\,\phi)=y_k(\phi)=\frac{\exp(a_k)}{\sum_j \exp(a_j)},
\hseq a_k=\trans{w_k}\phi
$$
で与えたときに最尤法を用いて直接$w_k$を求めよう.
$$
\dif{a_k}y_k=\frac{\exp(a_k)\left(\sum \exp(a_j)\right)-\exp(a_k)\exp(a_k)}{\left(\sum \exp(a_j)\right)^2}=y_k-y_k^2.
$$
$k\ne j$として
$$
\dif{a_j}y_k=-\frac{\exp(a_k) \exp(a_j)}{(\sum \exp(a_j))^2}=-y_k y_j.
$$
よってこの二つをまとめて
\begin{eqnarray}\label{ch4_mclass}
\dif{a_j}y_k=y_k(\delta_{kj}-y_j).
\end{eqnarray}

さて, 目的変数ベクトル$t_n$を$k$番目の要素だけが1であるものとし,
$(t_n)_k=t_{nk}$, $y_{nk}=y_k(\phi_n)$, $T=(t_{nk})$とおくと,
$$
p(T\,|\,w_1, \ldots, w_k)=\prod_{n=1}^N \prod_{k=1}^K p(C_k\,|\,\phi_n)^{t_{nk}}=\prod_{n,k} y_{nk}^{t_{nk}}.
$$
交差エントロピー誤差関数は
$$
E=-\log p(T\,|\,w_1, \ldots, w_k)=-\sum_{n,k} t_{nk} \log y_{nk}.
$$
よって
\begin{eqnarray*}
\dif{w_j}E
 &=& -\sum_{n,k} t_{nk} \frac{y_k(\phi_n)(\delta_{kj}-y_j(\phi_n))}{y_k(\phi_n)}\phi_n
 = -\sum_{n,k} t_{nk}(\delta_{kj} -y_{nj})\phi_n
\\
 &=& - \sum_n \left(\left(\sum_k t_{nk} \delta_{kj}\right)-\left(\sum_k t_{nk}\right)y_{nj}\right)\phi_n
 = -\sum_n (t_{nj}-y_{nj})\phi_n
\\
 &=& \sum(y_{nj}-t_{nj})\phi_n.
\end{eqnarray*}
やはり誤差$y_{nj}-t_{nj}$と基底関数$\phi_n$の積となる.

ヘッシアンをみる.
$$
\dif{w_k}y_{nj}=\dif{w_k}y_j(\phi_n)=y_k(\phi_n)(\delta_{kj}-y_j(\phi_n))\phi_n=y_{nk}(\delta_{kj}-y_{nj})\phi_n
$$
より
$$
H = \ddif{w_k}{w_j}E
  = \sum_n y_{nk}(\delta_{kj}-y_{nj})\phi_n \trans{\phi_n}.
$$

$H$が正定値であることを示そう.
任意の$M\times K$次元ベクトルを$u=\trans{(\trans{u_1}, \ldots, \trans{u_K})}$, $u_k$は$M$次元ベクトルとする.
$v_{nk}=\trans{u_k}\phi_n$, $f(x)=x^2$とおく. $f(x)$は下に凸.
\begin{eqnarray*}
\trans{u}H u
 &=& \sum_{n,k,j} y_{nk}(\delta_{kj}-y_{nj})(\trans{u_k}\phi_n)(\trans{\phi_n}u_j)\\
 &=& \sum_n\left(\sum_{k,j} y_{nk}\delta_{kj} v_{nk} v_{nj} - \sum_{k,j} y_{nk} y_{nj} v_{nk} v_{nj}\right)\\
 && （{\text 一般に}\left(\sum_k x_k\right)^2=\left(\sum_k x_k\right)\left(\sum_j x_j\right) = \sum_{k,j} x_k x_j{\text より}）\\
 &=& \sum_n \left(\sum_k y_{nk}v_{nk}^2 - \left(\sum_k y_{nk} v_{nk}\right)^2\right).
\end{eqnarray*}
ここで$\sum_k y_{nk}=1$, $0 < y_{nk} < 1$よりJensenの不等式を適用すると
$\trans{u}H u \ge 0$.

\section{プロビット回帰}
指数型分布族で表される条件付き確率分布に対して, クラスの事後確率はある線形関数とロジスティック（またはソフトマックス）関数の合成で表された.
$a=\trans{w}\phi$, $f(a)$を活性化関数として
$$
p(t=1\,|\,a)=f(a).
$$
と書ける範囲でもう少し考察する. $f(a)$がある確率密度$p(\theta)$の累積分布関数で表されるとする.
とくに$p(\theta)=\calN(\theta\,|\,0,1)$のとき累積分布関数は
$$
\Phi(a)=\int_{-\infty}^a \calN(\theta\,|\,0,1)\,d\theta.
$$
この逆関数をプロビット関数（probit）という. 誤差関数を
$$
\makeop{erf}(a)=\frac{2}{\sqrt{\pi}}\int_0^a \exp\left(-\theta^2\right) d\theta
$$
で定義する.
$x=\theta/\sqrt{2}$とおくと$dx=d\theta/\sqrt{2}$.
\begin{eqnarray*}
\Phi(a)
 &=& \int_{-\infty}^a = \int_{-\infty}^0 + \int_0^a=\half+\int_0^a \frac{1}{\sqrt{2\pi}}\exp\left(-\frac{\theta^2}{2}\right) d\theta\\
 &=&\half+\int_0^{a/\sqrt{2}}\frac{1}{\sqrt{2\pi}}\exp\left(-x^2\right)\sqrt{2}\,dx\\
 &=&\half\left(1+\frac{2}{\sqrt{\pi}}\int_0^{a/\sqrt{2}} \exp\left(-x^2\right) dx\right)
 =\half\left(1+\makeop{erf}\left(\frac{a}{\sqrt{2}}\right)\right).
\end{eqnarray*}
プロビット活性化関数を用いた一般化線形モデルをプロビット回帰という.

$x\rightarrow\infty$でロジスティック回帰の微分は$\sigma'(x)=\exp(-x)/(1+\exp(-x))^2 \sim \exp(-x)$.
プロビット回帰の微分は$\Phi'(x) \sim \exp(-x^2)$. つまりプロビット関数の逆関数はロジスティック関数よりも急速に1に近づき平らになる. プロビット回帰の方が外れ値に敏感.

\section{正準連結関数}
活性化関数として正準連結関数（canonical link function）と呼ばれるものを使い, 条件付き確率分布に指数型分布族を選んだときに誤差関数の微分が
「誤差」$\times$「特徴ベクトル」という形で書けることを示そう.
$$
p(t\,|\,\eta,s)=\frac{1}{s}h(t/s)g(\eta)\exp\left(\frac{\eta t}{s}\right)
$$
とする. 確率なので$\int p(t\,|\,\eta,s)\,dt=1$.
つまり
$$
g(\eta)\int h(t/s) \exp\left(\frac{\eta t}{s}\right) dt=s.
$$
$\eta$で微分して
\begin{eqnarray*}
&& \left(\dif{\eta}g(\eta)\right) \int h(t/s)\exp\left(\frac{\eta t}{s}\right) dt+g(\eta)\int (t/s) h(t/s) \exp\left(\frac{\eta t}{s}\right) dt\\
&&= \left(\dif{\eta}g(\eta)\right)\left(\frac{s}{g(\eta)}\right)+\int t p(t\,|\,\eta,s)\,dt
= s\dif{\eta}\log g(\eta) + \EE[t] = 0.
\end{eqnarray*}
よって
$$
y=\EE[t]=-s\dif{\eta} \log g(\eta).
$$
$y$が$\eta$の関数として表せた.
この逆関数が存在するとしてそれを$\eta=\psi(y)$と書くことにする.
$y$を連結関数$f(a)$と$w$の線形関数の合成,
$$
y=f(\trans{w}\phi)
$$
と書けるモデルを考える. 対数尤度関数は
$$
\log p(t\,|\,\eta,s)=\sum_{n=1}^N \log p(t_n\,|\,\eta,s)=\sum_{n=1}^N \left(\log g(\eta_n) + \frac{\eta_n t_n}{s}\right) + \text{const.}
$$
を考える. ここで$s$と$\eta$は独立, $\eta_n=\phi(y_n)$, $y_n=f(a_n)$, $a_n=\trans{w}\phi_n$.

パラメータが多いので依存関係に注意して微分する.
\begin{eqnarray*}&&
\dif{w}\eta_n=\psi'(y_n)f'(a_n)\phi_n,
\\&&
\dif{w}\log g(\eta_n)=\frac{g'(\eta_n)}{g(\eta_n)}\psi'(y_n)f'(a_n)\phi_n=-\frac{y_n}{s}\phi'(y_n)f'(a_n)\phi_n.
\end{eqnarray*}
よって
$$
\dif{w} \log p(t\,|\,\eta,s)=\sum_n \frac{1}{s}(t_n-y_n)\psi'(y_n)f'(a_n)\phi_n.
$$
連結関数として$f^{-1}(y)=\psi(y)$となるものを使ってみよう. $f(\psi(y))=y$を$y$で微分して
$$
f'(\psi(y))\psi'(y)=1.
$$
同じことだが$a=f^{-1}(y)=\psi(y)$を使って
$$
f'(a)\psi'(y)=1.
$$
よって
$$
\dif{w}E(w)=-\dif{w} \log p(t\,|\,\eta,s)=\frac{1}{s}\sum_n (y_n - t_n)\phi_n.
$$
「誤差」$\times$「特徴ベクトル」という形で書けることが分かった.

\section{ラプラス近似}
ロジスティック回帰のベイズ的な扱いは解析的に難しい.
ここではラプラス近似というものを紹介する. これは連続変数上の確率密度分布をあるガウス分布で近似することである(当然複数の山があると辛い).
$$
p(z)=\frac{1}{Z}f(z)
$$
$Z=\int f(z)\,dz$は正規化係数で未知とする. $p(z)$のモード, つまり最大値を与える$z_0$を探そう.
暗に山形を仮定しているので$f(z_0)>0$, $f''(z)<0$とする.
$f'(z_0)=0$だから$\log f(z)$を$z=z_0$の付近でテイラー展開すると
\begin{eqnarray*}
\log f(z) &\approx& \log f(z_0) + \frac{f'(z_0)}{f(z_0)}(z-z_0)+\half\frac{d^2}{dz^2} \log f(z_0) (z-z_0)^2\\
 &=& \log f(z_0) + \half\frac{d^2}{dz^2} \log f(z_0) (z-z_0)^2.
\end{eqnarray*}
$$
A=\left.-\frac{d^2}{dz^2}\log f(z) \right|_{z=z_0}
$$
とおくと
$$
\log f(z) \approx \log f(z) - \half A(z-z_0)^2.
$$
つまり
$$
f(z) \approx f(z_0) \exp\left(-\half A(z-z_0)^2\right).
$$
よって正規分布で近似すると
$$
q(z)=\left(\frac{A}{2\pi}\right)^{1/2} \exp\left(-\half A(z-z_0)^2\right).
$$

$M$次元の場合を考える. $p(z)=(1/Z)f(z)$. $z=z_0$が最大値を与えるなら$p(z_0)>0$, $\dif{z}f(z_0)=0$.
1次元と同様に
$$
A=\left.-\ddif{z_i}{z_j} \log f(z) \right|_{z=z_0}
$$
とすると$A>0$（正定値）で
$$
f(z) \approx f(z_0) \exp\left(-\half\trans{(z-z_0)}A(z-z_0)\right).
$$
よって
$$
q(z)=\calN(z\,|\,z_0,A^{-1})=\sqrt{\frac{|A|}{(2\pi)^M}} \exp\left(-\half\trans{(z-z_0)}A(z-z_0)\right).
$$
山が複数ある多峰的なときはどのモードを選ぶかでラプラス近似は異なる.
総的的にデータ数が多くなるとガウス分布に近づくので近似はよくなるが, ある点での近傍の情報しか利用していないため大域的な特徴がとらえられるとは限らない.

\section{モデルの比較とBIC}\label{ch4_bic}
前節のラプラス近似を行うと正規化係数$Z$の近似も分かる. ガウス分布の特性から
\begin{eqnarray}\label{ch4_approx}
Z&=&\int f(z)\,dx
 \approx \int f(z_0) \exp\left(-\half\trans{(z-z_0)}A(z-z_0)\right) dx
 = f(z_0)\sqrt{\frac{(2\pi)^M}{|A|}}.
 \qquad
\end{eqnarray}

データ集合$D$, パラメータ$\{\theta_i\}$の集合$\{M_i\}$を考えて各モデルに対して$p(D\,|\,\theta_i, M_i)$を定義する.
事前確率$p(\theta_i\,|\,M_i)$を決めてモデルエビデンス$p(D\,|\,M_i)$を計算してみよう. 以下$M_i$を略す.
また行列作用素$(\ddif{\theta_i}{\theta_j})$を$\nabla^2$と書くことにする.
$$
p(D)=\int p(D\,|\,\theta)\,p(\theta)\,d\theta.
$$
$f(\theta)=p(D\,|\,\theta)\,p(\theta)$, $Z=p(D)$とすると$\theta=\theta_{\text{MAP}}$のときのラプラス近似を用いて
$$
A=-\nabla^2 \log p(D\,|\,\theta_{\text{MAP}})\,p(\theta_{\text{MAP}}).
$$
\begin{eqnarray*}
\log p(D)&=&\log Z \approx \log f(\theta_{\text{MAP}}) + \log \sqrt{\frac{(2\pi)^M}{|A|}}\\
         &=&\log p(D\,|\,\theta_{\text{MAP}})+\left(\log p(\theta_{\text{MAP}})+\frac{M}{2} \log (2\pi) - \half \log |A|\right).
\end{eqnarray*}
括弧内の後ろ三項をOccam係数という.

この値をごくごく粗く近似してみる.
$$
p(\theta)=\calN(\theta\,|\,m,B^{-1}),
$$
$$
H=-\nabla^2\log p(D\,|\,\theta_{\text{MAP}})
$$
とすると
$$
A=H-\nabla^2\log p(\theta_{\text{MAP}})=H+B.
$$
$$
\log p(\theta)=\log|B|^{1/2}-(M/2)\log(2\pi)-(1/2)\trans{(\theta-m)}B{(\theta-m)}
$$
を使って
$$
\log p(D) \approx \log p(D\,|\,\theta_{\text{MAP}})-\half\trans{(\theta_{\text{MAP}}-m)}B(\theta_{\text{MAP}}-m)-\half\log |H+B|+\half\log |B|.
$$
データ集合$D$の点の個数を$N$とし, それぞれが独立であると考えると$H$の各要素は$O(N)$. $B$は$N$や次元の数$M$には依存しないので$M$, $N$を大きくしたとき無視できるとすると$C$を定数として
$$
\log |H| \approx \log \prod^M (N C) = M \log N + M \log C \approx M \log N.
$$
よって
$$
\log p(D) \approx \log p(D\,|\,\theta_{\text{MAP}}) - \half M \log N.
$$
これをベイズ情報量基準（Bayesian Information Criterion, BIC）と呼ばれるモデルの良さを評価するための指標に一致する.
かなり無理筋な近似ではあるが, ラプラス近似が（別の方法で導出される）BICと関連があることを暗示している.

\section{ディラックのデルタ関数}
ヘヴィサイド関数$H(x)$を$\RR$から$\RR$への関数で
$$
H(x)=
\begin{cases}
1 & (x>0), \\
0 & (x \le 0)
\end{cases}
$$
と定義する. 気持ち的にはデルタ関数はヘヴィサイド関数の微分である.
$$
H'(x) = \delta(x).
$$
$x=0$以外では微分すると$0$で, $x=0$では無限大になるので
$$
\delta(x)=
\begin{cases}
\infty & (x=0), \\
0 & (x \neq 0)
\end{cases}
$$
という感じである.
ただ, 実際にはデルタ関数は積分を通してしか扱われない. 厳密には測度論や関数解析の理論を用いて正当化されなければならないが
次のような関係式が成り立つ.
$$
\int_{-\infty}^\infty \delta(x)\,dx = \Bigl[H(x)\Bigr]_{-\infty}^\infty = H(\infty) - H(-\infty) = 1.
$$
実数値関数$f(x)$に対して部分積分を適用すると
\begin{eqnarray*}
\int_{-\infty}^\infty f(x) \delta(x)\,dx &=& \Bigl[H(x)f(x)\Bigr]_{-\infty}^\infty - \int_{-\infty}^\infty H(x) f'(x)\,dx\\
 &=& f(\infty) - \int_0^\infty f'(x)\,dx
 = f(\infty) - (f(\infty) - f(0))
 = f(0).
\end{eqnarray*}
\vspace{0pt}

\section{ロジスティックシグモイド関数とプロビット関数の逆関数}
\begin{itemize}
\item[] ロジスティックシグモイド
$$\sigma(a) = \frac{1}{1+e^{-a}}.$$
\item[] プロビット関数の逆関数
$$\Phi(a) = \int_{-\infty}^a \calN(\theta\,|\,0,1)\,d\theta.$$
\end{itemize}
$\Phi(\lambda a)$で$\sigma(a)$を近似するように$\lambda$を調節する.
近似にあたっては, これらが原点で同じ傾きをもつようにしよう.
\pagebreak
\begin{eqnarray*}&&
\sigma'(0)
=\sigma(a)(1-\sigma(a))\Bigr|_{a=0}
=\half\left(1-\half\right)=\frac{1}{4},
\\&&
\left.\dif{a}\left(\Phi'(\lambda a)\right)\right|_{a=0}
=\left.\frac{1}{\sqrt{2\pi}}\lambda
 \exp\left(-\half(\lambda a)^2\right)\right|_{a=0}
=\frac{\lambda}{\sqrt{2\pi}}.
\end{eqnarray*}
この二つが等しいので
$
\sfrac{\lambda}{\sqrt{2\pi}}=\sfrac{1}{4}
$
より$\lambda=\sqrt{\pi/8}$を得る.

$\Phi$と$\calN$に関する畳み込み計算の関係式：$\lambda>0$として
\begin{equation}\label{inv_probit}
\int_{-\infty}^\infty \Phi(\lambda a)\,\calN(a\,|\,\mu,\sigma^2)\,da=\Phi\left(\frac{\lambda \mu}{\sqrt{1+\lambda^2\sigma^2}}\right).
\end{equation}
これを示そう. 左辺を$L$, 右辺を$R$とおく.
まずガウス分布に関する積分を簡単にする.
$$
\calN(a\,|\,\mu,\sigma^2)=\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{1}{2\sigma^2}(x-\mu)^2\right).
$$
$$
f(x)=\calN(x\,|\,0,1)=\frac{1}{2\pi}\exp\left(-\half x^2\right)
$$
とおく.
$$
\Phi(x)=\int_{-\infty}^x f(y)\,dy.
$$
$a=\mu+\sigma x$とおくと$da=\sigma dx$で
$$
\calN(a\,|\,\mu,\sigma^2)=\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{1}{2\sigma^2}(\mu+\sigma x-\mu)^2\right)=\frac{1}{\sigma}\frac{1}{\sqrt{2\pi}}\exp\left(-\half x^2\right)=\frac{f(x)}{\sigma}.
$$
よって
$$
L=\int_{-\infty}^\infty \Phi(\lambda \mu+\lambda\sigma x)\frac{f(x)}{\sigma}\sigma\,dx=\int_{-\infty}^\infty\Phi(\lambda \mu + \lambda \sigma x)f(x)\,dx.
$$
まず$L$と$R$の$\mu$に関する微分が等しいことを示す. $\Phi'(x)=f(x)$なので
\begin{eqnarray*}
\dif{\mu}R&=&f\left(\frac{\lambda \mu}{\sqrt{1+\lambda^2\sigma^2}}\right) \frac{\lambda}{\sqrt{1+\lambda^2\sigma^2}},
\\
\dif{\mu}L&=&\int_{-\infty}^\infty \lambda f(\lambda \mu + \lambda \sigma x)f(x)\,dx
\\
&=& \lambda \int_{-\infty}^\infty \left(\frac{1}{\sqrt{2\pi}}\right)^2\exp\left(-\half(\lambda\mu+\lambda \sigma x)^2-\half x^2\right) dx.
\end{eqnarray*}
$\exp()$の内側を$x$について平方完成しよう：
\begin{eqnarray*}
&&-\half\left((\lambda\sigma)^2+1\right)x^2+2\lambda^2\mu \sigma x + \lambda^2 \mu^2)\\
&&\quad=-\half\left(\left(\sqrt{1+\lambda^2\sigma^2}x+\frac{\lambda^2\mu \sigma}{\sqrt{1+\lambda^2\sigma^2}}\right)^2+\lambda^2 \mu^2-\frac{\lambda^4 \mu^2 \sigma^2}{1+\lambda^2 \sigma^2} \right).
\end{eqnarray*}
定数項は
$$
\lambda^2\mu^2-\frac{\lambda^4 \mu^2 \sigma^2}{1+\lambda^2 \sigma^2}=\frac{\lambda^2\mu^2+\lambda^4\mu^2\sigma^2-\lambda^4\mu^2\sigma^2}{1+\lambda^2 \sigma^2}=\frac{\lambda^2\mu^2}{1+\lambda^2 \sigma^2},
$$
$$
\int_{-\infty}^\infty\exp\left(-\frac{1}{2\sigma^2}x^2\right) dx=\sqrt{2\pi}\sigma
$$
より
$$
I:=\int_{-\infty}^\infty \exp\left(-\half\left(\sqrt{1+\lambda^2\sigma^2}x+\frac{\lambda^2\mu\sigma}{\sqrt{1+\lambda^2\sigma^2}}\right)^2\right) dx=\frac{\sqrt{2\pi}}{\sqrt{1+\lambda^2\sigma^2}}.
$$
よって
\begin{eqnarray*}
\dif{\mu}L
 &=& \frac{\lambda}{2\pi} \exp\left(-\half\left(\frac{\lambda^2\mu^2}{1+\lambda^2 \sigma^2}\right)\right) I
 = \frac{\lambda}{2\pi} \exp\left(-\half\left(\frac{\lambda^2\mu^2}{1+\lambda^2 \sigma^2}\right)\right) \frac{\sqrt{2\pi}}{\sqrt{1+\lambda^2\sigma^2}}\\
 &=& \frac{\lambda}{\sqrt{1+\lambda^2\sigma^2}}f\left(\frac{\lambda\mu}{\sqrt{1+\lambda^2\sigma^2}}\right)
 = \dif{\mu}R.
\end{eqnarray*}
つまり$L$と$R$は定数の差を除いて等しいことが分かった. 定数項が$0$であることを示す. $\mu \rightarrow \infty$で
\begin{eqnarray*}
&& L \rightarrow \int_{-\infty}^\infty \Phi(\infty)f(x)\,dx=\int_{-\infty}^\infty f(x)\,dx=1,
\\
&& R \rightarrow \int_{-\infty}^\infty f(y)\,dy=1.
\end{eqnarray*}
よって$L=R$が示された.

\section{ベイズロジスティック回帰}\label{ch4_bayes}
ベイズロジスティック回帰にラプラス近似を行ってみよう.

事前確率分布：$y_n=\sigma(\trans{w}\phi_n)$とおいて
$$
p(t\,|\,w)=\prod_{n=1}^N y_n^{t_n}(1-y_n)^{1-t_n}.
$$
事前ガウス分布：ハイパーパラメータ$m_0$, $S_0$を用いて
$$
p(w)=\calN (w\,|\,m_0, S_0).
$$
事後分布は$N$次元縦ベクトル$t=\trans{(t_0, \ldots, t_N)}$を用いて
$$
p(w\,|\,t) \propto p(w)\,p(t\,|\,w).
$$
よって
$$
\log p(w\,|\,t)
= -\half\trans{(w-m_0)}S_0^{-1}(w-m_0)
  +\sum_{n=1}^N \left(t_n \log y_n {+} (1-t_n)\log (1-y_n)\right) + \text{const.}
$$
これを最大化する最大事後確率（MAP）の解$w_\text{MAP}$を何らかの方法で求める.
ラプラス近似をするために$w_\text{MAP}$を平均とするガウス分布で近似しよう. 以後$w_\text{MAP}$を$m_N$と書く.
共分散は今まで何度もやった通り
$$
S_N^{-1}=-\nabla^2 \log p(w\,|\,t) = S_0^{-1} + \sum_{n=1}^N y_n(1-y_n)\phi_n \trans{\phi_n}.
$$
よって事後確率分布をガウス分布で近似すると
$$
q(w)=\calN(w\,|\,m_N, S_N).
$$
次に$C_1$についての予測分布を求める.
$$
p(C_1\,|\,\phi,t)=\int p(C_1\,|\,\phi,w)\,p(w\,|\,t)\,dw \approx \int \sigma(\trans{w}\phi)\,q(w)\,dw.
$$
デルタ関数の性質から
$$
\sigma(\trans{w}\phi)=\int \delta(a-\trans{w}\phi)\,\sigma(a)\,da.
$$
よって
\begin{eqnarray*}
p(C_1\,|\,\phi,t)
 &\approx& \int \left(\int \delta(a-\trans{w}\phi)\,\sigma(a)\,q(w)\,da\right) dw\\
 &=&\int \left(\int \delta(a-\trans{w}\phi)\,q(w)\,dw\right) \sigma(a)\,da
 =\int p(a)\,\sigma(a)\,da.
\end{eqnarray*}
ただし
$$
p(a)=\int \delta(a-\trans{w}\phi)\,q(w)\,dw
$$
とおいた.
$p(a)$の平均を求める. $q(w)$の平均が$m_N$であることに注意すると,
\begin{eqnarray*}
\mu_a
 &=& \EE[a]
  =\int p(a)a\,da
  =\int \int \delta(a-\trans{w}\phi)\,q(w)a\,dw da\\
 &=& \int \left(\int \delta(a-\trans{w}\phi)a\, da\right) q(w)\,dw
  = \int q(w)(\trans{w}\phi)\,dw\\
 &=& \trans{\left(\int q(w)w\,dw\right)}\phi = \trans{\EE[w]}\phi = \trans{m_N}\phi.
\end{eqnarray*}
分散は
\begin{eqnarray*}
\sigma(a)^2 &=& \int p(a)(a^2-\EE[a]^2)\,da \\
 &=& \int (\delta(a-\trans{w}\phi)a^2\,da)\,q(w)\,dw
     - \int (\delta(a-\trans{w}\phi)\EE[a]^2\,da)\,q(w)\,dw\\
 &=&\int q(w)(\trans{w}\phi)^2\,dw - \int q(w) (\trans{m_N}\phi)^2\,dw
 \\
 &=& \trans{\phi}\left( \int q(w)(w\trans{w}-m_N \trans{m_N})\,dw \right)\phi.
\end{eqnarray*}
括弧内は
$$
\EE[w\trans{w}]-m_N\trans{m_N}\int q(w)\,dw=(m_N\trans{m_N}+S_N)-m_N \trans{m_N}=S_N
$$
より
$$
\sigma_a^2=\trans{\phi}S_N\phi.
$$
これはPRML式(3.59)に示されている線形回帰モデルの予測分布の分散
$$
\sigma_N^2(x)=\frac{1}{\beta}+\trans{\phi(x)}S_N\phi(x)
$$
で, $\beta \to \infty$としてノイズを消したものは$\sigma_a^2$に一致する.
よって予測分布の近似は
$$
p(C_1\,|\,t)=\int \sigma(a)\,p(a)\,da=\int \sigma(a) \,\calN(a\,|\,\mu_a,\sigma_a^2)\,da.
$$

別の方法でも求めてみる：
$w$の座標を座標変換して第一成分の単位基底ベクトル$e$が$\phi$と同じ向きになるようにとり, 残りは$\phi$と直交するようにとる.
つまり
$$
e=\phi/\Norm{\phi}.
$$
それに関する$w$の係数を$w=\trans{(w_1,w_2)}$と書く.
$w_1$は1次元ベクトルの係数で非負, $w_2$は$M-1$次元である.
すると$w$と$\phi$の内積は$e$の方向のみが残るので
$$
\trans{w}\phi=w_1\Norm{\phi}.
$$
$q(w)=q(w_1,w_2)=q(w_2\,|\,w_1)\,q(w_1)$より
\begin{eqnarray*}
p &=& \int \sigma(\trans{w}\phi)\,q(w)\,dw=\int \sigma(w_1\Norm{\phi}) \int q(w_2\,|\,w_1)\,dw_2 dw_1\\
 &=& \int \sigma(w_1\Norm{\phi}) \left( \int q(w_2\,|\,w_1)\,dw_2 \right) q(w_1)\,dw_1\\
 && \text{（括弧内は$1$なので）}\\
 &=& \int \sigma(w_1\Norm{\phi})\,q(w_1)\,dw_1.
\end{eqnarray*}
$q(w)\,\calN(w\,|\,m_N,S_N)$で$w=\trans{(w_1,w_2)}$と書いたときの
\pagebreak
$w_1$に関する周辺分布は$e$方向への射影化になる：
$$
q(w_1)=\calN(w_1\,|\,\trans{e}m_N,\trans{e}S_N e).
$$
$a=w_1\Norm{\phi}$と書くと$a$に関しては平均は$\Norm{\phi}$倍, 分散は$\Norm{\phi}^2$倍されるので
$$
q(a)
= \fnvPair
  \calN(a|\trans{(\Norm{\phi}e)}m_N,\trans{\Norm{\phi}e}S_N \Norm{\phi}e)
= \fnvPair
  \calN(a|\trans{\phi}m_N,\trans{\phi}S_N \phi).
$$
ロジスティックシグモイド関数でのガウス分布の畳み込み積分は解析的には難しい.
ここではロジスティックシグモイド関数をプロビット関数の逆関数で近似することで畳み込み積分を解析的に扱おう.
前節の結果から$\lambda = \sqrt{\pi/8}$として$\sigma(a) \approx \Phi(\lambda a)$と近似すると
\begin{eqnarray}\label{ch4_probit}
\int \sigma(a)\,\calN(a\,|\,\mu,\sigma^2)\,da &\approx& \sigma\left(\frac{\mu}{\sqrt{1+\lambda^2 \sigma^2}}\right) \nonumber \\
&& （\kappa(\sigma^2)=1/\sqrt{1+(\pi/8)\sigma^2}{\text とおくと}）\nonumber \\
&=&\sigma\left(\kappa(\sigma^2)\mu\right).
\end{eqnarray}
よって
$$
p(C_1\,|\,\phi,t) \approx \sigma(\kappa(\sigma_a^2)\,\mu_a)
$$
という近似予測分布が得られた.