はじめに
PRML解答例まとめを参照
演習 2.41
ガンマ関数
\Gamma(a) \equiv \int_{0}^{\infty} u^{a-1} e^{-u} \mathrm{d} u \tag{1.141}
の定義から,ガンマ分布
\operatorname{Gam}(\lambda | a, b)=\frac{1}{\Gamma(a)} b^{a} \lambda^{a-1} \exp (-b \lambda) \tag{2.146}
が正規化されていることを示せ.
ガンマ分布が正規化されている\displaystyle \left(\int_0^{\infty}\frac{1}{\Gamma(a)} b^{a} \lambda^{a-1} \exp (-b \lambda) d\lambda = 1\right)ことを示す。
これにはb\lambda = uという変数変換を使うと(d\lambda = b^{-1}du)、
\begin{aligned}
&\int_0^{\infty}\frac{1}{\Gamma(a)} b^{a} \lambda^{a-1} \exp (-b \lambda) d\lambda \\
=&\frac{1}{\Gamma(a)}\int_0^{\infty}b^a\lambda^{a-1}\exp(-u)b^{-1}du \\
=&\frac{1}{\Gamma(a)}\int_0^{\infty}u^{a-1}e^{-u}du = \frac{\Gamma(a)}{\Gamma(a)} = 1
\end{aligned}
となるので、正規化されていることが示された。
演習 2.42
ガンマ分布
\operatorname{Gam}(\lambda | a, b)=\frac{1}{\Gamma(a)} b^{a} \lambda^{a-1} \exp (-b \lambda) \tag{2.146}
の平均,分散,およびモードを求めよ.
まず,平均を求める.
\begin{aligned}
\mathbb E \left[\lambda \right] =& \int_0^{\infty} \frac{1}{\Gamma(a)} b^{a} \lambda^{a-1} \exp (-b \lambda) \lambda \mathrm{d} \lambda \\
=& \frac{1}{\Gamma(a)} \int_0^{\infty} b^{a} \lambda^a \exp (-b \lambda) \mathrm{d} \lambda \\
\end{aligned}
ここでb\lambda = uと変数変換すると,\mathrm{d} \lambda = b^{-1} \mathrm{d} uであるから
\begin{aligned}
\mathbb E \left[\lambda \right] =& \frac{1}{\Gamma(a)} \int_0^{\infty} u^a \exp (-u) b^{-1} \mathrm{d} u \\
=& \frac{1}{\Gamma(a)} b^{-1} \int_0^{\infty} u^a e^{-u} \mathrm{d} u \\
=& \frac{1}{\Gamma(a)} b^{-1} \Gamma(a+1) \\
=& \frac{1}{\Gamma(a)} b^{-1} a\Gamma(a) \\
=& \frac{a}{b}
\end{aligned}
となる.
次に,分散を求める.
\begin{aligned}
\mathbb E \left[\lambda^2 \right] =& \int_0^{\infty} \frac{1}{\Gamma(a)} b^{a} \lambda^{a-1} \exp (-b \lambda) \lambda^2 \mathrm{d} \lambda \\
=& \frac{1}{\Gamma(a)} \int_0^{\infty} b^{a} \lambda^{a+1} \exp (-b \lambda) \mathrm{d} \lambda \\
\end{aligned}
ここでb\lambda = uと変数変換すると,\mathrm{d} \lambda = b^{-1} \mathrm{d} uであるから
\begin{aligned}
\mathbb E \left[\lambda^2 \right] =& \frac{1}{\Gamma(a)} \int_0^{\infty} b^{-1} u^{a+1} \exp (-u) b^{-1} \mathrm{d} u \\
=& \frac{1}{\Gamma(a)} b^{-2} \int_0^{\infty} u^{a+1} e^{-u} \mathrm{d} u \\
=& \frac{1}{\Gamma(a)} b^{-2} \Gamma(a+2) \\
=& \frac{1}{\Gamma(a)} b^{-2} (a+1)\Gamma(a+1) \\
=& \frac{1}{\Gamma(a)} b^{-2} a(a+1)\Gamma(a) \\
=& \frac{a(a+1)}{b^2}
\end{aligned}
となる.したがって分散は
\begin{aligned}
\operatorname{var}\left[\lambda \right] =& \mathbb E \left[\lambda^2 \right] - \mathbb E \left[\lambda \right]^2 \\
=& \frac{a(a+1)}{b^2} - \frac{a^2}{b^2} \\
=& \frac{a}{b^2}
\end{aligned}
となる.
次に,モードを求める.モードは分布の停留点を求めれば良いので(2.146)を\lambdaで微分して0とおく.
\begin{aligned}
& \frac{\mathrm{d}}{\mathrm{d}\lambda} \operatorname{Gam}(\lambda | a, b) = 0 \\
\Leftrightarrow ~ & \frac{\mathrm{d}}{\mathrm{d}\lambda} \frac{1}{\Gamma(a)} b^{a} \lambda^{a-1} \exp (-b \lambda) = 0 \\
\Leftrightarrow ~ & \frac{\mathrm{d}}{\mathrm{d}\lambda} \lambda^{a-1} e^{-b \lambda} = 0 \\
\Leftrightarrow ~ & (a-1) \lambda^{a-2} e^{-b \lambda} + \lambda^{a-1} e^{-b \lambda} (-b) = 0 \\
\Leftrightarrow ~ & \lambda^{a-2} e^{-b \lambda} (a-1-b\lambda) = 0 \\
\Leftrightarrow ~ & \lambda = \frac{a-1}{b} \\
\end{aligned}
したがってモードは\displaystyle \frac{a-1}{b}となる.
演習 2.43
次の分布は,1変数ガウス分布を一般化したものである.
p\left(x | \sigma^{2}, q\right)=\frac{q}{2\left(2 \sigma^{2}\right)^{1 / q} \Gamma(1 / q)} \exp \left(-\frac{|x|^{q}}{2 \sigma^{2}}\right) \tag{2.293}
この分布が,次のように正規化されていることを示せ.
\int_{-\infty}^{\infty} p\left(x | \sigma^{2}, q\right) \mathrm{d} x=1 \tag{2.294}
そして,q=2で(2.293)がガウス分布となることを示せ.目的変数がt=y(\mathbf{x}, \mathbf{w})+\epsilonで,ランダムノイズ変数\epsilonが分布(2.293)に従う回帰モデルを考える.入力ベクトル集合\mathbf{X} = \{ \mathbf{x}_1, \cdots, \mathbf{x}_N \}と,目的変数に相当する観測値が\mathbf{t} = (t_1, \ldots, t_N)^{\mathrm{T}}であるとき,\mathbf{w}と\sigma^2についての対数尤度関数が
\ln p\left(\mathbf{t} | \mathbf{X}, \mathbf{w}, \sigma^{2}\right)=-\frac{1}{2 \sigma^{2}} \sum_{n=1}^{N}\left|y\left(\mathbf{x}_{n}, \mathbf{w}\right)-t_{n}\right|^{q}-\frac{N}{q} \ln \left(2 \sigma^{2}\right)+\mathrm{const} \tag{2.295}
になることを示せ.ただし,\mathrm{const}は\mathbf{w}と\sigma^2の両方と独立な項である.なお,これは\mathbf{w}の関数として見た場合,1.5.5節で扱ったL_q誤差関数である.
まず(2.293)について
\begin{aligned}
\int_{-\infty}^{\infty} p(x | \sigma^2, q)dx &= \int_{-\infty}^{\infty} \frac{q}{2\left(2 \sigma^{2}\right)^{1 / q} \Gamma(1 / q)} \exp \left(-\frac{|x|^{q}}{2 \sigma^{2}}\right) d x
\\
&=\frac{q}{2\left(2 \sigma^{2}\right)^{1 / q} \Gamma(1 / q)} \cdot 2 \int_{0}^{\infty} \exp \left(-\frac{|x|^{q}}{2 \sigma^{2}}\right) d x \\
&=\frac{q}{\left(2 \sigma^{2}\right)^{1 / q} \Gamma(1 / q)} \int_{0}^{\infty} \exp \left(-\frac{x^{q}}{2 \sigma^{2}}\right) d x
\end{aligned}
ここで\displaystyle u = \frac{x^q}{2\sigma^2}とおくと、\displaystyle du = \frac{q}{2\sigma^2}x^{q-1}dx, x=(2\sigma^2u)^{1/q}となるので、
\begin{aligned}
\int_{0}^{\infty} \exp \left(-\frac{x^{2}}{2 \sigma^{2}}\right) d x
&=\int_{0}^{\infty} \frac{2 \sigma^{2}}{q x^{q-1}} e^{-u} du \\
&=\frac{2 \sigma^{2}}{q} \int_{0}^{\infty} \frac{1}{\left(2 \sigma^{2} u\right)^{\frac{q-1}{q}}} e^{-u} d u \\
&=\frac{2 \sigma^{2}}{q} \cdot \frac{\left(2 \sigma^{2}\right)^{\frac{1}{q}}}{2 \sigma^{2}} \int_{0}^{\infty} u^{\frac{1}{q}-1} e^{-u} d u \\
&=\frac{\left(2 \sigma^{2}\right)^{\frac{1}{q}}}{q} \Gamma(1 / q)
\end{aligned}
よって
\int_{-\infty}^{\infty} p\left(x \mid \sigma^{2}, q\right) d x=\frac{q}{\left(2 \sigma^{2}\right)^{1 / q} \Gamma(1 / q)} \cdot \frac{\left(2 \sigma^{2}\right)^{\frac{1}{q}}}{q} \Gamma(1 / q)=1
正規化されていることが示された。
また、q=2のとき\Gamma(1/2) = \sqrt{\pi}を代入すれば
\int_{-\infty}^{\infty} p\left(x \mid \sigma^{2}, 2 \right) dx = \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left( -\frac{x^2}{2\sigma^2}\right)
となり、ガウス分布の形になる。
最後に、ランダムノイズが\epsilonがp(x|\sigma^2,q)に従うとすると、\epsilon = t-y(\mathbf{x}, \mathbf{w})が(2.293)式のxに相当する。
\mathbf{X} = \{\mathbf{x}_1, \ldots, \mathbf{x}_N\}を入力として目的変数に相当する観測値が\mathbf{t} = (t_1, \ldots, t_N)^{\mathrm{T}}とするとき、1.2.5節での議論から
p(\mathbf{t} | \mathbf{X}, \mathbf{w}, \sigma^2) = \prod_{n=1}^N p(t_n | \mathbf{x}_n, \mathbf{w}, \sigma^2)
となるので、これの対数尤度は
\begin{aligned}
\ln p\left(\mathbf{t} | \mathbf{X}, \mathbf{w}, \sigma^2 \right) &=\ln \left\{\prod_{n=1}^{N} \frac{q}{2\left(2 \sigma^{2}\right)^{1 / q} \Gamma(1 / q)} \exp \left(-\frac{\mid t_n-y\left(\mathbf{x}_n, \mathbf{w} \right) \mid^q}{2 \sigma^{2}}\right)\right \} \\
&=-\frac{1}{2 \sigma^{2}} \sum_{n=1}^{N}\left|y\left(\mathbf{x}_n, \mathbf{w} \right)-t_{n}\right|^{q}-\frac{N}{q} \ln \left(2 \sigma^{2}\right)+\operatorname{const}
\end{aligned}
となる(\textrm{const}は\mathbf{w}と\sigma^2の両方に独立な項である)。
演習 2.44
1変数ガウス分布\mathcal{N}\left(x | \mu, \tau^{-1}\right)について考える.共役事前分布はガウス–ガンマ分布
p(\mu, \lambda)=\mathcal{N}\left(\mu \mid \mu_{0},(\beta \lambda)^{-1}\right) \operatorname{Gam}(\lambda \mid a, b) \tag{2.154}
で,独立同分布な観測値集合が\mathbf{x}=\left\{x_{1}, \ldots, x_{N}\right\}であるとする.事後分布も,事前分布と同じガウスーガンマ分布になることを示し,各パラメータに対する事後分布の式を書き下せ.
事後分布はp(\mu,\lambda|\mathbf{X})、尤度関数はp(\mathbf{X}|\mu, \lambda)と書ける。
このとき共役事前分布p(\mu,\lambda)と合わせて、ベイズの定理から
p(\mu,\lambda|\mathbf{X}) \propto p(\mathbf{X}|\mu, \lambda)p(\mu,\lambda)
と書ける。また、(2.152)式のように尤度関数は
\begin{aligned}
p(\mathbf{X} \mid \mu, \lambda)&=\prod_{n=1}^{N}\left(\frac{\lambda}{2 \pi}\right)^{1 / 2} \exp \left\{-\frac{\lambda}{2}\left(x_{n}-\mu\right)^{2}\right\} \\
&\propto\left[\lambda^{1 / 2} \exp \left(-\frac{\lambda \mu^{2}}{2}\right)\right]^{N} \exp \left\{\lambda \mu \sum_{n=1}^{N} x_{n}-\frac{\lambda}{2} \sum_{n=1}^{N} x_{n}^{2}\right\}
\end{aligned}
の形で書ける。これと
p(\mu, \lambda) \propto (\beta\lambda)^{1/2}\exp\left( -\frac{\beta\lambda}{2}(\mu-\mu_0)^2 \right)\lambda^{a-1}\exp(-b\lambda)
を利用して、事後分布p(\mu,\lambda|\mathbf{X})の指数部を計算してみてガウス-ガンマ分布のそれと同型になることを示せば良い(正規化定数を求めるのはキツイ)。
\begin{aligned}
p(\mu,\lambda|\mathbf{X}) &= \lambda^{\frac{N}{2}}(\beta \lambda)^{\frac{1}{2}} \lambda^{a-1} \exp \left\{-\frac{N \lambda \mu^{2}}{2}-\frac{\beta \lambda}{2}\left(\mu-\mu_{0}\right)^{2}-b \lambda+\lambda \mu \sum_{n=1}^{N} x_{n}-\frac{\lambda}{2} \sum_{n=1}^{N} x_{n}^{2}\right\} \\
&= \beta^{\frac{1}{2}} \lambda^{\frac{N}{2}+\frac{1}{2}+a-1} \exp \left\{\mu^{2} \left( -\frac{N \lambda}{2}-\frac{\beta \lambda}{2}\right)+\mu\left(\beta \lambda \mu_{0}+\lambda \sum_{n=1}^{N} x_{n}\right)-b \lambda-\frac{\lambda}{2} \sum_{n=1}^{N} x_{n}^{2}-\frac{\beta \lambda}{2} \mu_{0}^{2}\right\} \\
&=\beta^{\frac{1}{2}} \lambda^{\frac{N}{2}+\frac{1}{2}+a-1} \exp \left\{ -\frac{\lambda(N+\beta)}{2} \mu^{2}+\mu \lambda\left(\beta \mu_{0}+\sum_{n=1}^{N} x_{n}\right)-\lambda\left(b+\frac{1}{2} \sum_{n=1}^{N} x_{n}^{2}-\frac{\beta}{2} \mu_{0}^{2}\right)\right\} \\
&=\beta^{\frac{1}{2}} \lambda^{\frac{N}{2}+\frac{1}{2}+a-1} \exp \left\{-\frac{\lambda(N+\beta)}{2}\left(\mu-\frac{\beta\mu_0+\sum_{n=1}^{N} x_{n}}{N+\beta}\right)^{2}\right\}
\exp \left\{-\lambda\left(b+\frac{1}{2} \sum_{n=1}^N x_{n}^{2}+\frac{\beta}{2} \mu_{0}^{2}-\frac{\left(\beta\mu_0+\sum_{n=1}^{N} x_{n}\right)^{2}}{2(N+\beta)}\right)\right\} \\
&\propto (\lambda(N+\beta))^{\frac{1}{2}}\exp \left\{-\frac{\lambda(N+\beta)}{2}\left(\mu-\frac{\beta\mu_0+\sum_{n=1}^{N} x_{n}}{N+\beta}\right)^{2}\right\}\lambda^{a+\frac{N}{2}-1}\exp \left\{-\left(b+\frac{1}{2} \sum_{n=1}^N x_{n}^{2}+\frac{\beta}{2} \mu_{0}^{2}-\frac{\left(\beta\mu_0+\sum_{n=1}^{N} x_{n}\right)^{2}}{2(N+\beta)}\right)\lambda\right\}
\end{aligned}
これは
\begin{aligned}
\mu_{N} &=\frac{\beta \mu_{0}+\sum_{n=1}^{N} x_{n}}{N+\beta} \\
\lambda_{N}^{-1} &= (\lambda(N+\beta))^{-1} \\
a_{N} &=a+\frac{N}{2} \\
b_{N} &=b+\frac{1}{2} \sum_{n=1}^N x_{n}^{2}+\frac{\beta}{2} \mu_{0}^{2}-\frac{\left(\beta\mu_0+\sum_{n=1}^{N} x_{n}\right)^{2}}{2(N+\beta)}
\end{aligned}
としたときの\mathcal{N}(\mu|\mu_N,\lambda_N^{-1})\operatorname{Gam}(\lambda|a_N, b_N)の指数部分と同型になる。
したがって、事後分布もガウス-ガンマ分布の形になっていることが示された。
演習 2.45
\mathcal{W}(\mathbf{\Lambda} \mid \mathbf{W}, \nu)=B|\mathbf{\Lambda}|^{(\nu-D-1) / 2} \exp \left(-\frac{1}{2} \operatorname{Tr}\left(\mathbf{W}^{-1} \mathbf{\Lambda}\right)\right) \tag{2.155}
で定義されたウィシャート分布が確かに,多変量ガウス分布の精度行列の共役事前分布であることを確かめよ.
事前分布と尤度と掛け合わせて事後分布を求めたとき,目的とする変数についての形式が変わらないような事前分布が共役事前分布である.したがって,尤度を求め,精度\mathbf \Lambdaについてウィシャート分布の形になっていることを確認すれば良い.
\mathbf{X} = \{\mathbf{x}_1, \mathbf{x}_2, \cdots,\mathbf{x}_N\}が与えられたとき,精度\mathbf \Lambdaのガウス分布\mathcal{N}(\mathbf{X}|\boldsymbol{\mu},\mathbf{\Lambda}^{-1})の尤度関数は
\begin{aligned}
& \prod_{n=1}^N \mathcal{N}(\mathbf{x}_n|\boldsymbol{\mu},\mathbf{\Lambda}^{-1}) \\
\propto & |\mathbf{\Lambda}|^{\frac{1}{2}} \exp{\left( -\frac{1}{2} \sum_{n=1}^N (\mathbf{x}_n - \boldsymbol{\mu})^{\mathrm{T}} \mathbf{\Lambda} (\mathbf{x}_n - \boldsymbol{\mu}) \right)} \\
=& |\mathbf{\Lambda}|^{\frac{1}{2}} \exp{\left( -\frac{1}{2} \sum_{n=1}^N \operatorname{Tr} \left[(\mathbf{x}_n - \boldsymbol{\mu})^{\mathrm{T}} \mathbf{\Lambda} (\mathbf{x}_n - \boldsymbol{\mu}) \right]\right)} \\
=& |\mathbf{\Lambda}|^{\frac{1}{2}} \exp{\left( -\frac{1}{2} \sum_{n=1}^N \operatorname{Tr} \left[(\mathbf{x}_n - \boldsymbol{\mu}) (\mathbf{x}_n - \boldsymbol{\mu})^{\mathrm{T}} \mathbf{\Lambda} \right]\right)} \\
=& |\mathbf{\Lambda}|^{\frac{1}{2}} \exp{\left( -\frac{1}{2} \operatorname{Tr} \left[ \sum_{n=1}^N (\mathbf{x}_n - \boldsymbol{\mu}) (\mathbf{x}_n - \boldsymbol{\mu})^{\mathrm{T}} \mathbf{\Lambda} \right]\right)} \\
=& |\mathbf{\Lambda}|^{\frac{1}{2}} \exp{\left( -\frac{1}{2} \operatorname{Tr} \left[\mathrm{S} \mathbf{\Lambda} \right]\right)} \\
\end{aligned}
となり,ウィシャート分布と同じ形であることがわかる.ただし\displaystyle \sum_{n=1}^N (\mathbf{x}_n - \boldsymbol{\mu}) (\mathbf{x}_n - \boldsymbol{\mu})^{\mathrm{T}} = \mathrm{S}とおいた.
以上より,ウィシャート分布(2.156)は多変量ガウス分布の精度行列の共役事前分布となっていることが確認できた.
演習 2.46
\begin{align}
p(x|\mu,a,b)&=\int_{0}^{\infty} \mathcal{N} (x | \mu, \tau^{-1}) \operatorname{Gam}(\tau | a, b) d \tau \\
&=\int_{0}^{\infty} \frac{b^{a} e^{(-b \tau)} \tau^{a-1}}{\Gamma(a)}\left(\frac{\tau}{2 \pi}\right)^{1 / 2} \exp \left\{-\frac{\tau}{2}(x-\mu)^{2}\right\} \mathrm{d} \tau\\
&=\frac{b^{a}}{\Gamma(a)}\left(\frac{1}{2 \pi}\right)^{\frac{1}{2}} \int_{0}^{\infty} \tau^{a-\frac{1}{2}} \exp \left\{-\tau\left(b+\frac{(x-\mu)^{2}}{2}\right)\right\} d \tau \tag{2.146}
\end{align}
の積分を計算し
\operatorname{St}(x \mid \mu, \lambda, \nu)=\frac{\Gamma(\nu / 2+1 / 2)}{\Gamma(\nu / 2)}\left(\frac{\lambda}{\pi \nu}\right)^{1 / 2}\left[1+\frac{\lambda(x-\mu)^{2}}{\nu}\right]^{-\nu / 2-1 / 2} \tag{2.159}
になることを確かめよ.
まず(2.158)式に\displaystyle \operatorname{Gam} ( \lambda |a, b)=\frac{1}{r(a)} b^{a} \lambda^{a-1} \exp (-b \lambda), \quad \Gamma(x)=\int_{0}^{\infty} u^{x-1} e^{-u} dxを適用する。
\begin{aligned}
&\int_{0}^{\infty} \mathcal{N} (x | \mu, \tau^{-1}) \operatorname{Gam}(\tau | a, b) d \tau \\
=&\int_{0}^{\infty} \frac{\sqrt{\tau}}{\sqrt{2 \pi}} \exp \left\{ -\frac{\tau}{2}(x-\mu)^{2}\right\} \frac{b^{a}}{\Gamma(a)} \tau^{a-1} \exp (-b \tau) d \tau \\
=&\frac{b^{a}}{\Gamma(a)}\left(\frac{1}{2 \pi}\right)^{\frac{1}{2}} \int_{0}^{\infty} \tau^{a-\frac{1}{2}} \exp \left\{-\tau\left(b+\frac{(x-\mu)^{2}}{2}\right)\right\} d \tau
\end{aligned}
ここで、\displaystyle b+\frac{(x-\mu)^2}{2}=cとおき、\tau c = uとすると、d\tau = c^{-1}du
\begin{aligned}
=&\frac{b^{a}}{\Gamma(a)}\left(\frac{1}{2 \pi}\right)^{\frac{1}{2}} \int_{0}^{\infty}\left(\frac{u}{c}\right)^{a-\frac{1}{2}} \exp (-u) \cdot c^{-1} d u \\
=&\frac{b^{a}}{\Gamma(a)}\left(\frac{1}{2 \pi}\right)^{\frac{1}{2}} c^{-a-\frac{1}{2}} \int_{0}^{\infty} u^{a-\frac{1}{2}} e^{-u} d u \\
=&\frac{b^{a}}{\Gamma(a)}\left(\frac{1}{2 \pi}\right)^{\frac{1}{2}}\left[b+\frac{(x-\mu)^{2}}{2}\right]^{-a-\frac{1}{2}} \Gamma\left(a+\frac{1}{2}\right)
\end{aligned}
よって(2.158)の式変形が示された。次に\displaystyle \nu=2a,\quad \lambda=\frac{a}{b}とおくと, \displaystyle b^a = \left( \frac{\nu}{2\lambda} \right)^\frac{\nu}{2}= \left( \frac{2\lambda}{\nu} \right)^{- \frac{\nu}{2}}
\begin{aligned}
&=\frac{\Gamma\left(\frac{\nu}{2}+\frac{1}{2}\right)}{\Gamma\left(\frac{\nu}{2}\right)}\left(\frac{1}{2 \pi}\right)^{\frac{1}{2}}\left( \frac{2\lambda}{\nu} \right)^{- \frac{\nu}{2}}\left[\frac{\nu}{2 \lambda}+\frac{(x-\mu)^{2}}{2}\right]^{-\frac{\nu}{2}-\frac{1}{2}} \\
&=\frac{\Gamma\left(\frac{\nu}{2}+\frac{1}{2}\right)}{\Gamma\left(\frac{\nu}{2}\right)} \left(\frac{\lambda}{\pi \nu}\right)^{\frac{1}{2}} \left(\frac{2 \lambda}{\nu} \cdot \frac{\nu+\lambda(x-\mu)^{2}}{2 \lambda}\right)^{-\frac{\nu}{2}-\frac{1}{2}} \\
&=\frac{\Gamma\left(\frac{\nu}{2}+\frac{1}{2}\right)}{\Gamma\left(\frac{\nu}{2}\right)}\left(\frac{\lambda}{\pi \nu}\right)^{\frac{1}{2}}\left[1+\frac{\lambda(x-\mu)^{2}}{\nu}\right]^{-\frac{\nu}{2}-\frac{1}{2}}
\end{aligned}
となり、(2.159)式に変形できることが示された。
演習 2.47
\nu \to \inftyの極限で,t分布
\operatorname{St}(x \mid \mu, \lambda, \nu)=\frac{\Gamma(\nu / 2+1 / 2)}{\Gamma(\nu / 2)}\left(\frac{\lambda}{\pi \nu}\right)^{1 / 2}\left[1+\frac{\lambda(x-\mu)^{2}}{\nu}\right]^{-\nu / 2-1 / 2} \tag{2.159}
がガウス分布になることを示せ.ヒント: 正規化係数を無視し,xへの依存性だけに注目する.
スチューデントのt分布が\nu \to \inftyの極限で\displaystyle \mathcal{N}(x|\mu, \lambda^{-1})=\sqrt{\frac{\lambda}{2\pi}} \exp \left\{ -\frac{\lambda(x-\mu)^2}{2} \right\}になることを示す。
\displaystyle \left[1+\frac{\lambda(x-\mu)^{2}}{\nu}\right]^{-\frac{\nu}{2}-\frac{1}{2}}の部分を\exp(\cdot)で表すように変形させることを意識する。また、\displaystyle \frac{\Gamma\left(\frac{\nu}{2}+\frac{1}{2}\right)}{\Gamma\left(\frac{\nu}{2}\right)\left(\frac{\nu}{2}\right)^{1/2}}が\nu \to \inftyの極限で1に収束することを示す。このために\displaystyle \lim_{n\to\infty}\frac{\Gamma(x+n)}{\Gamma(n)n^x} = 1であることを示す(※公式の解答集や問題文のヒントから、ここはどうせ規格化定数で吸収するので考えなくても良いことになっているけれど、一応示す)。
まずガンマ関数の定義から
\Gamma(x) = \lim_{n\to\infty}\frac{(n-1)!n^x}{\prod_{k=0}^{n-1}(x+k)}
逆数をとって
\frac{1}{\Gamma(x)}=\lim_{n\to\infty}\frac{\prod_{k=0}^{n-1}(x+k)}{(n-1)!n^x}
よって
\begin{aligned}
1 &= \lim_{n\to\infty}\Gamma(x)\frac{\prod_{k=0}^{n-1}(x+k)}{(n-1)!n^x} \\
&= \lim_{n\to\infty}\frac{\Gamma(x)\Gamma(x+n)}{\Gamma(n)\Gamma(x)n^x} \\
&= \lim_{n\to\infty}\frac{\Gamma(x+n)}{\Gamma(n)n^x}
\end{aligned}
ここで、x = 1/2, n=\nu/2とすれば、\displaystyle \lim_{n\to\infty}\frac{\Gamma\left(\frac{\nu}{2}+\frac{1}{2}\right)}{\Gamma\left(\frac{\nu}{2}\right)\left(\frac{\nu}{2}\right)^{1/2}} = 1となることが示された。
ちなみにガンマ関数の定義をWikipediaにならって
\Gamma(x) = \lim_{n\to\infty}\frac{n!n^x}{\prod_{k=0}^{n}(x+k)}
とした場合、上の式変形を進めていくと
1 = \lim_{n\to\infty}\left( 1+\frac{x}{n+1}\right)\frac{\Gamma(x+n)}{\Gamma(n)n^x}
となる。n\to\inftyで\displaystyle 1+\frac{x}{n+1} \to 1となることを使えば結果は同様となる(ということであってるかな……?)。
続いて\displaystyle \left[1+\frac{\lambda(x-\mu)^{2}}{\nu}\right]^{-\frac{\nu}{2}-\frac{1}{2}} について
\left[1+\frac{\lambda(x-\mu)^{2}}{\nu}\right]^{-\frac{\nu}{2}-\frac{1}{2}}= \exp \left\{ -\frac{\nu + 1}{2}\ln \left( 1+ \frac{\lambda(x-\mu)^2}{\nu}\right)\right\}
ここでテイラー展開\ln (1+\epsilon) = \epsilon+O(\epsilon^2)を用いると
\begin{aligned}
&\exp \left\{ -\frac{\nu + 1}{2}\ln \left( 1+ \frac{\lambda(x-\mu)^2}{\nu}\right)\right\} \\
=&\exp\left\{ -\frac{\nu+1}{\nu}\cdot \left( \frac{\lambda(x-\mu)^2}{2}+O(\nu^{-2})\right)\right\} \\
=&\exp\left\{ -\frac{\nu+1}{\nu}\cdot \frac{\lambda(x-\mu)^2}{2}+O(\nu^{-1})\right\}
\end{aligned}
\nu \to \inftyの極限で、上式は\displaystyle \exp\left\{ -\frac{\lambda(x-\mu)^2}{2}\right\}となる。
したがって以上をまとめると、\nu \to \inftyの極限で
\begin{aligned}
&\lim_{\nu\to\infty}\frac{\Gamma\left(\frac{\nu}{2}+\frac{1}{2}\right)}{\Gamma\left(\frac{\nu}{2}\right)}\left(\frac{\lambda}{\pi \nu}\right)^{\frac{1}{2}}\left[1+\frac{\lambda(x-\mu)^{2}}{\nu}\right]^{-\frac{\nu}{2}-\frac{1}{2}} \\
=&\lim_{\nu\to\infty}\frac{\Gamma\left(\frac{\nu}{2}+\frac{1}{2}\right)}{\Gamma\left(\frac{\nu}{2}\right)\left(\frac{\nu}{2}\right)^\frac{1}{2}} \left(\frac{\lambda}{2\pi}\right)^{\frac{1}{2}}\left[1+\frac{\lambda(x-\mu)^{2}}{\nu}\right]^{-\frac{\nu}{2}-\frac{1}{2}} \\
=&1\cdot\left(\frac{\lambda}{2\pi}\right)^{\frac{1}{2}}\exp\left\{ -\frac{\lambda(x-\mu)^2}{2}\right\} \\
=&\mathcal{N}(x|\mu,\lambda^{-1})
\end{aligned}
となる。よってガウス分布となることが示された。
別解
ネイピア数の定義(\lim_{n\to\infty}(1 + n^{-1})^{n}=e)に基づいて示す方法
\begin{aligned}
\lim_{\nu\to\infty}\left[1+\frac{\lambda(x-\mu)^{2}}{\nu}\right]^{-\frac{\nu+1}{2}}=& \lim_{\nu\to\infty}\left\{\left[1+\frac{\lambda(x-\mu)^{2}}{\nu}\right]^{\frac{\nu}{\lambda(x-\mu)^{2}}}\right\}^{\frac{\lambda(x-\mu)^{2}}{\nu}({-\frac{\nu+1}{2}})}\\
=&\displaystyle \exp\left\{ -\frac{\lambda(x-\mu)^2}{2}\right\}
\end{aligned}
{} の中身は \nu\to\infty で exp になり {} の外側の指数部分は \nu\to\infty で \frac{\lambda(x-\mu)^{2}}{2} となることを用いた
演習 2.48
1変数のスチューデントのt分布
\operatorname{St}(x \mid \mu, \lambda, \nu)=\frac{\Gamma(\nu / 2+1 / 2)}{\Gamma(\nu / 2)}\left(\frac{\lambda}{\pi \nu}\right)^{1 / 2}\left[1+\frac{\lambda(x-\mu)^{2}}{\nu}\right]^{-\nu / 2-1 / 2} \tag{2.159}
の導出で用いたのと同様の手続きで,スチューデントのt分布の多変量形式である
\operatorname{St}(\mathbf{x} \mid \boldsymbol{\mu}, \mathbf{\Lambda}, \nu)=\frac{\Gamma(D / 2+\nu / 2)}{\Gamma(\nu / 2)} \frac{|\mathbf{\Lambda}|^{1 / 2}}{(\pi \nu)^{D / 2}}\left[1+\frac{\Delta^{2}}{\nu}\right]^{-D / 2-\nu / 2} \tag{2.162}
の結果を,
\operatorname{St}(\mathbf{x} \mid \boldsymbol{\mu}, \mathbf{\Lambda}, \nu)=\int_{0}^{\infty} \mathcal{N}\left(\mathbf{x} \mid \boldsymbol{\mu},(\eta \mathbf{\Lambda})^{-1}\right) \operatorname{Gam}(\eta \mid \nu / 2, \nu / 2) \mathrm{d} \eta \tag{2.161}
の変数\etaを周辺化することで確かめよ.定義(2.161)を用いて,積分変数の順序を交換することで,多変量t分布が正規化されていることを示せ.
(2.161)の積分表現から(2.162)を得ることを目的とする。途中で以下の定理を用いる。
任意のn\times n行列\mathbf{A}と任意のスカラー値kに対して
|k\mathbf{A}| = k^n|\mathbf{A}|
が成り立つ(統計のための行列代数P.217, 系13.2.4)
途中でマハラノビス距離\Delta^2 = (\mathbf{x}-\boldsymbol{\mu})^{\mathrm{T}}\mathbf{\Lambda}(\mathbf{x}-\boldsymbol{\mu})を用いると
\begin{aligned}
\operatorname{St}(\mathbf{x} \mid \boldsymbol{\mu}, \mathbf{\Lambda}, v)
&=\int_{0}^{\infty} \mathcal{N}\left(\mathbf{x} \mid \boldsymbol{\mu},(\eta \mathbf{\Lambda})^{-1}\right) \operatorname{Gam}\left(\eta \left| \frac{\nu}{2}, \frac{\nu}{2}\right.\right) d \eta \\
&=\int_{0}^{\infty} \frac{1}{(2 \pi)^{D/2}}|\eta \mathbf{\Lambda}|^{\frac{1}{2}} \exp \left\{-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^{\mathrm{T}}(\eta \mathbf{\Lambda})(\mathbf{x}-\boldsymbol{\mu})\right\} \frac{1}{\Gamma\left(\frac{\nu}{2}\right)}\left(\frac{\nu}{2}\right)^{\frac{\nu}{2}} \eta^{\frac{\nu}{2}-1} \exp \left\{-\frac{\nu}{2} \eta\right\} d\eta \\
&=\frac{|\mathbf{\Lambda}|^{\frac{1}{2}}}{(2 \pi)^{D/2}} \frac{\left(\frac{\nu}{2}\right)^{\nu/2}}{\Gamma\left(\frac{\nu}{2}\right)} \int_{0}^{\infty} \eta^{\frac{D}{2}} \cdot \eta^{\frac{\nu}{2}-1} \exp \left\{-\frac{\nu}{2} \eta-\frac{\eta}{2} \Delta^{2}\right\} d\eta
\end{aligned}
ここで\displaystyle \tau = \left( \frac{\nu}{2}+\frac{\Delta^2}{2} \right)とおくと、\displaystyle d\tau = \left( \frac{\nu}{2}+\frac{\Delta^2}{2}\right)d\etaで、
\begin{aligned}
(与式)
&=\frac{|\mathbf{\Lambda}|^{\frac{1}{2}}}{(2 \pi)^{D/2}} \frac{\left(\frac{\nu}{2}\right)^{\nu/2}}{\Gamma\left(\frac{\nu}{2}\right)}\left(\frac{\nu}{2}+\frac{\Delta^{2}}{2}\right)^{-\frac{D}{2}-\frac{\nu}{2}} \int_{0}^{\infty} \tau^{\frac{D}{2}+\frac{\nu}{2}-1} e^{-\tau} d \tau \\
&=\frac{|\mathbf{\Lambda}|^{\frac{1}{2}}}{(2 \pi)^{D/2}} \frac{\left(\frac{\nu}{2}\right)^{\nu/2}}{\Gamma\left(\frac{\nu}{2}\right)}\left(\frac{\nu}{2}+\frac{\Delta^{2}}{2}\right)^{-\frac{D}{2}-\frac{\nu}{2}} \Gamma\left(\frac{D}{2}+\frac{\nu}{2}\right) \\
&=\frac{\Gamma\left(\frac{D}{2}+\frac{\nu}{2}\right)}{\Gamma\left(\frac{\nu}{2}\right)} \cdot \frac{\left(\frac{\nu}{2}\right)^{\nu/2}}{(2 \pi)^{D/2}} | \mathbf{\Lambda} |^{\frac{1}{2}}\left(\frac{\nu}{2}\right)^{-\frac{D}{2}-\frac{\nu}{2}}\left[1+\frac{\Delta^{2}}{\nu}\right]^{-\frac{D}{2}-\frac{\nu}{2}} \\
&=\frac{\Gamma\left(\frac{D}{2}+\frac{\nu}{2}\right)}{\Gamma\left(\frac{\nu}{2}\right)} \frac{|\mathbf{\Lambda}|^{\frac{1}{2}}}{(\nu \pi)^{D/2}}\left[1+\frac{\Delta^{2}}{\nu}\right]^{-\frac{D}{2}-\frac{\nu}{2}}
\end{aligned}
よって(2.162)式が得られた。
また、この多変数t分布が正規化されていることを示す。すなわち、
\int_{-\infty}^{\infty}\operatorname{St}(\mathbf{x}|\boldsymbol{\mu},\mathbf{\Lambda},\nu) d\mathbf{x} = 1
であることを示す。(問題文のヒントから)積分が交換可能であることを利用すると、定義から
\begin{aligned}
& \int_{-\infty}^{\infty} \operatorname{St}(\mathbf{x}|\boldsymbol{\mu},\mathbf{\Lambda},\nu) d\mathbf{x} \\
=& \int_{-\infty}^{\infty} \int_{0}^{\infty} \mathcal{N}\left(\mathbf{x} | \boldsymbol{\mu},\left(\eta\mathbf{\Lambda}^{-1}\right)\right) \operatorname{Gam}\left(\eta \left| \frac{\nu}{2}, \frac{\nu}{2}\right.\right) d\eta d\mathbf{x} \\
=& \int_{0}^{\infty} \underbrace{\int_{-\infty}^{\infty} \mathcal{N}\left(\mathbf{x} | \boldsymbol{\mu},\left(\eta \mathbf{\Lambda}^{-1}\right)\right) d\mathbf{x}}_{1} \operatorname{Gam}\left(n \left| \frac{\nu}{2}, \frac{\nu}{2}\right.\right) d \eta \\
=& \int_{0}^{\infty} \operatorname{Gam}\left(\eta \left| \frac{\nu}{2}, \frac{\nu}{2}\right.\right) d\eta \\
=& 1
\end{aligned}
ガンマ分布が正規化されていることは演習問題2.41で示した。
演習 2.49
ガウス分布とガンマ分布のたたみ込みである多変量スチューデントt分布の定義(2.161)を用いて,(2.162)で定義される多変量t分布の平均,共分散,およびモード
\mathbb{E}[\mathbf{x}] = \boldsymbol{\mu} \hspace{2em} (\mu > 1のとき)\tag{2.164}
\operatorname{cov}[\mathbf{x}] = \frac{\nu}{\nu -2} \mathbf{\Lambda}^{-1} \hspace{2em} (\mu > 2のとき)\tag{2.165}
\operatorname{mode}[\mathbf{x}] = \boldsymbol{\mu} \tag{2.166}
を確かめよ.
\begin{aligned}
\mathbb{E}[\mathbf{x}] &= \int_{0}^{\infty} \underbrace{\int
\mathcal{N}\left(\mathbf{x} | \boldsymbol{\mu},\left(\eta\mathbf{\Lambda}\right)^{-1}\right) \mathbf{x} d\mathbf{x}}_{\boldsymbol{\mu}(\because{(2.58)式})}\operatorname{Gam}\left(\eta \left| \frac{\nu}{2}, \frac{\nu}{2}\right.\right) d\eta \\
&= \int_{0}^{\infty} \boldsymbol{\mu}\operatorname{Gam}\left(\eta \left| \frac{\nu}{2}, \frac{\nu}{2}\right.\right) d\eta \\
&=\boldsymbol{\mu} (\because 演習2.41)
\end{aligned}
\begin{aligned}
\operatorname{cov}[\mathbf{x}] &= \int\int_{0}^{\infty}
\mathcal{N}\left(\mathbf{x} | \boldsymbol{\mu},\left(\eta\mathbf{\Lambda}\right)^{-1}\right) \operatorname{Gam}\left(\eta \left| \frac{\nu}{2}, \frac{\nu}{2}\right.\right) d\eta (\mathbf{x}-\boldsymbol{\mu})(\mathbf{x}-\boldsymbol{\mu})^\mathrm{T} d\mathbf{x} \\
&= \int_{0}^{\infty}\underbrace{\int
\mathcal{N}\left(\mathbf{x} | \boldsymbol{\mu},\left(\eta\mathbf{\Lambda}\right)^{-1}\right) (\mathbf{x}-\boldsymbol{\mu})(\mathbf{x}-\boldsymbol{\mu})^\mathrm{T} d\mathbf{x}}_{\left(\eta\mathbf{\Lambda}\right)^{-1}} \operatorname{Gam}\left(\eta \left| \frac{\nu}{2}, \frac{\nu}{2}\right.\right) d\eta \\
&=\int_{0}^{\infty}\left(\eta\mathbf{\Lambda}\right)^{-1} \frac{1}{\Gamma(\frac{\nu}{2})}\left(\frac{\nu}{2}\right)^\frac{\nu}{2}\eta^{\left(\frac{\nu}{2}-1\right)} e^{{-\frac{\nu}{2}}\eta} d\eta \\
&=\mathbf{\Lambda}^{-1} \frac{\frac{\nu}{2}^\frac{\nu}{2}}{\Gamma(\frac{\nu}{2})} \int_{0}^{\infty}\eta^{\left(\frac{\nu}{2}-2\right)}e^{-\frac{\nu}{2}\eta}d\eta \\\\
&ここで\frac{\nu}{2}\eta = uとおく。d\eta=\frac{2}{\nu}du \\\\
&=\mathbf{\Lambda}^{-1} \frac{\frac{\nu}{2}^\frac{\nu}{2}}{\Gamma(\frac{\nu}{2})} \int_{0}^{\infty}\left(\frac{2u}{\nu}\right)^{\left(\frac{\nu}{2}-2\right)}e^{-u}\frac{2}{\nu}du \\
&=\mathbf{\Lambda}^{-1} \frac{\frac{\nu}{2}^\frac{\nu}{2}}{\Gamma(\frac{\nu}{2})}\left(\frac{2}{\nu}\right)^{\left(\frac{\nu}{2}-1\right)} \underbrace{\int_{0}^{\infty}u^{\left(\frac{\nu}{2}-2\right)}e^{-u}du}_{\Gamma(\frac{\nu}{2}-1)} \\
&=\mathbf{\Lambda}^{-1} \frac{\Gamma(\frac{\nu}{2}-1)}{\Gamma(\frac{\nu}{2})}\frac{\nu}{2} \\
&=\mathbf{\Lambda}^{-1} \frac{\Gamma(\frac{\nu}{2}-1)}{\Gamma(\frac{\nu}{2}-1)\left(\frac{\nu}{2}-1\right)}\frac{\nu}{2}\\
&=\frac{\nu}{\nu-2}\Lambda^{-1}\\
\end{aligned}
mode[x]=\boldsymbol{\mu} を示す
\begin{aligned}\\
\operatorname{St}(\mathbf{x} \mid \boldsymbol{\mu}, \mathbf{\Lambda}, \nu)=\frac{\Gamma(D / 2+\nu / 2)}{\Gamma(\nu / 2)} \frac{|\mathbf{\Lambda}|^{1 / 2}}{(\pi \nu)^{D / 2}}\left[1+\frac{\Delta^{2}}{\nu}\right]^{-D / 2-\nu / 2}\\
\end{aligned}
このようにかけるので
\begin{aligned}
&\frac{\partial}{\partial\mathbf{x}}\operatorname{St}(\mathbf{x} \mid \boldsymbol{\mu}, \mathbf{\Lambda}, \nu)\propto \frac{\partial}{\partial\mathbf{x}}\left[1+\frac{\Delta^{2}}{\nu}\right]^{-D / 2-\nu / 2}\\
&=\left(-\frac{D}{2}-\frac{\nu}{2}\right)\left[1+\frac{\Delta^2}{\nu}\right]^{-D/2-\nu/2-1}\frac{2}{\nu}\mathbf{\Lambda}(\mathbf{x}-\boldsymbol{\mu})\\
\end{aligned}
これが\mathbf{0}となるのは\mathbf{x}=\boldsymbol{\mu}のとき
演習 2.50
\nu \to \inftyの極限で,多変量スチューデントt分布
\operatorname{St}(\mathbf{x} \mid \boldsymbol{\mu}, \mathbf{\Lambda}, \nu)=\frac{\Gamma(D / 2+\nu / 2)}{\Gamma(\nu / 2)} \frac{|\mathbf{\Lambda}|^{1 / 2}}{(\pi \nu)^{D / 2}}\left[1+\frac{\Delta^{2}}{\nu}\right]^{-D / 2-\nu / 2} \tag{2.162}
が,平均が\boldsymbol{\mu}で精度が\mathbf{\Lambda}のガウス分布になることを示せ.
\lim_{\nu \to \infty} \operatorname{St}(\mathbf{x} \mid \boldsymbol{\mu}, \mathbf{\Lambda}, \nu)=\mathcal{N}\left(\mathbf{x} | \boldsymbol{\mu},\mathbf{\Lambda}\right)を示す。スチューデントのt分布の形は
\operatorname{St}(\mathbf{x} \mid \boldsymbol{\mu}, \mathbf{\Lambda}, \nu)=\frac{\Gamma(D / 2+\nu / 2)}{\Gamma(\nu / 2)} \frac{|\mathbf{\Lambda}|^{1 / 2}}{(\pi \nu)^{D / 2}}\left[1+\frac{\Delta^{2}}{\nu}\right]^{-D / 2-\nu / 2}
ここで\displaystyle \lim_{\nu \to \infty}\frac{\Gamma(D/2+\nu/2)}{\Gamma(\nu/2)\left(\nu/2\right)^{D/2}}=1であるので(演習2.47を参照)、
\begin{aligned}\\
\lim_{\nu\to\infty}\operatorname{St}(\mathbf{x} \mid \boldsymbol{\mu}, \mathbf{\Lambda}, \nu) &= \lim_{\nu\to\infty}\frac{|\mathbf{\Lambda}|^{1 / 2}}{(2\pi)^{D / 2}}\left[1+\frac{\Delta^{2}}{\nu}\right]^{-D / 2-\nu / 2} \\
&=\frac{|\mathbf{\Lambda}|^{1 / 2}}{(2\pi)^{D / 2}}\lim_{\nu\to\infty}\underbrace{\left[1+\frac{\Delta^{2}}{\nu}\right]^{-D / 2}}_{\to1}\left[1+\frac{\Delta^{2}}{\nu}\right]^{-\nu / 2}
\end{aligned}
ここで\displaystyle u=-\frac{\nu}{2}と置換し\displaystyle \exp(x)=\lim_{n\to\infty}\left[1+\frac{x}{n}\right]^{n}=\lim_{n\to-\infty}\left[1+\frac{x}{n}\right]^{n}を用いると
\begin{aligned}
\lim_{\nu\to\infty}\operatorname{St}(\mathbf{x} \mid \boldsymbol{\mu}, \mathbf{\Lambda}, \nu) &= \frac{|\mathbf{\Lambda}|^{1 / 2}}{(2\pi)^{D / 2}}\lim_{u\to-\infty}\left[1-\frac{\Delta^{2}}{2u}\right]^{u}\\
&=\frac{|\mathbf{\Lambda}|^{1 / 2}}{(2\pi)^{D / 2}}\exp\left(-\frac{\Delta^{2}}{2}\right)\\
&=\mathcal{N}\left(\mathbf{x} | \boldsymbol{\mu},\mathbf{\Lambda}\right)
\end{aligned}
演習 2.51
本章の周期変数の議論で用いた,いろいろな三角関数の公式は,次の関係を用いて容易に証明できる.
\exp(iA) = \cos A+i\sin A \tag{2.296}
ただしiは-1の平方根である.
\cos^2 A+\sin^2 A = 1 \tag{2.177}
の結果を
\exp(iA) \exp(-iA)=1 \tag{2.297}
から証明せよ.同様に
\cos (A-B)=\Re \exp \{i(A-B)\} \tag{2.298}
を用いて
\cos A \cos B+\sin A \sin B=\cos (A-B) \tag{2.178}
を証明せよ.ただし,\Reは実部を示す.最後に,\sin (A-B) = \Im \exp\{i(A-B)\}から,
\sin (A-B)=\sin A \cos B-\cos A \sin B \tag{2.183}
を証明せよ.ただし,\Imは虚部を示す.
(2.296)を(2.297)に代入すると
\begin{aligned}( \text{左辺} ) &= \exp(iA)\exp(-iA) \\&= (\cos A + i \sin iA)(\cos A - i \sin iA) \\&= \cos^2 A + \sin^2 A \end{aligned}
よって、(2.177)は示された。
(2.298)より
\begin{aligned} \cos (A-B) &= \Re \exp \{ i(A-B) \} \\&= \Re \exp(iA) \exp(-iB) \\&= \Re (\cos A + i \sin iA)(\cos B - i \sin iB) \\&= \cos A \cos B + \sin A \sin B \end{aligned}
また
\begin{aligned} \sin (A-B) &= \Im \exp \{ i(A-B) \} \\&= \Im \exp(iA) \exp(-iB) \\&= \Im (\cos A + i \sin iA)(\cos B - i \sin iB) \\&= \sin A \cos B - \cos A \sin B \end{aligned}
演習 2.52
フォン・ミーゼス分布
p\left(\theta \mid \theta_{0}, m\right)=\frac{1}{2 \pi I_{0}(m)} \exp \left\{m \cos \left(\theta-\theta_{0}\right)\right\} \tag{2.179}
は,mが大きいとき,モード\theta_0の周囲で鋭く尖る.\xi = m^{1/2}(\theta-\theta_0)と定義し,余弦関数のテイラー展開が
\cos \alpha = 1-\frac{\alpha^2}{2}+O(\alpha^4) \tag{2.299}
であるとして,m \to \inftyでフォン・ミーゼス分布がガウス分布になることを示せ.
\xi = m^{1/2}(\theta-\theta_0)より、\theta-\theta_0 = \xi m^{-1/2}。
これをフォン・ミーゼス分布(2.179)に代入すると、
\begin{aligned} p\left(\theta \mid \theta_{0}, m\right) &\propto \exp \left\{m \cos \left(\theta-\theta_{0}\right)\right\} \\&= \exp \left\{m \cos \left(\xi m^{-1/2}\right)\right\} \\&= \exp \left\{m \left( 1 - \frac{1}{2} \xi^2 m^{-1} + O(\xi^4 m^{-2}) \right) \right\} \\&\propto \exp \left( -\frac{\xi^2}{2} \right) \\&= \exp \left\{ - \frac{m(\theta - \theta_{0})^2}{2} \right\} \end{aligned}
これはガウス分布の形になっているので示された。
演習 2.53
三角関数の公式
\sin (A-B)=\sin A \cos B-\cos A \sin B \tag{2.183}
を用いて,\theta_0についての
\sum_{n=1}^{N} \sin \left(\theta_{n}-\theta_{0}\right)=0 \tag{2.182}
の解が
\theta_{0}^{\mathrm{ML}}=\tan ^{-1}\left\{\frac{\sum_{n} \sin \theta_{n}}{\sum_{n} \cos \theta_{n}}\right\} \tag{2.184}
となることを示せ.
(2.183)で示されている加法定理\sin(A-B) = \sin A \cos B- \cos A \sin Bを用いて(2.182)式を解いて
\sum_{n=1}^{N}\sin(\theta_n-\theta_0)=0 \\
\sum_{n=1}^{N}\left( \sin \theta_n \cos \theta_0 -\cos \theta_n \sin\theta_0 \right) =0 \\
sin\theta_0\sum_{n=1}^{N}\cos(\theta_n)=cos\theta_0\sum_{n=1}^{N}\sin(\theta_n)
\tan \theta_0 = \frac{\sum_{n}\sin\theta_n}{\sum_{n}\cos\theta_n}
なので、最尤推定量\theta_0^{\mathrm{ML}}は
\theta_0^{\mathrm{ML}} = \tan^{-1}\left\{ \frac{\sum_{n}\sin\theta_n}{\sum_{n}\cos\theta_n} \right\}
となり、(2.184)式を得ることができた。
演習 2.54
フォン・ミーゼス分布
p\left(\theta \mid \theta_{0}, m\right)=\frac{1}{2 \pi I_{0}(m)} \exp \left\{m \cos \left(\theta-\theta_{0}\right)\right\} \tag{2.179}
の1階と2階の導関数を求め,さらにm>0でI_0(m) > 0であることを用いて,分布は\theta = \theta_0で最大になり,\theta = \theta_0 + \pi\, (\mathrm{mod}\ 2\pi)で最小になることを示せ.
1階の導関数は
\frac{\partial}{\partial \theta} p\left(\theta \mid \theta_{0}, m\right)=\frac{1}{2 \pi I_{0}(m)} \exp \left(m \cos \left(\theta-\theta_{0}\right)\right)\left(-m \sin \left(\theta-\theta_{0}\right)\right)
2階の導関数は
\begin{aligned}
\frac{\partial^{2}}{\partial \theta^{2}} p\left(\theta \mid \theta_{0}, m\right) &=
\frac{1}{2 \pi I_{0}(m)} \exp \left(m \cos\left(\theta-\theta_{0}\right)\left(-m \sin \left(\theta-\theta_{0}\right)\right)\left(-m \sin \left(\theta-\theta_{0}\right)\right)\right. \\ &\left.\left.+\exp \left(m\cos \left(\theta-\theta_{0}\right)\right)( -m\cos \left(\theta-\theta_{0}\right)\right)\right\} \\ &=\frac{1}{2 \pi I_{0}(m)} \exp \left(m\cos \left(\theta-\theta_{0}\right)\right)\left(m^{2} \sin ^{2}\left(\theta-\theta_{0}\right)-m \cos \left(\theta-\theta_{0}\right)\right)
\end{aligned}
上の結果より\displaystyle \frac{\partial p}{\partial \theta} = 0となるのは\displaystyle \frac{1}{2\pi I_0(m)}\exp (m \cos (\theta - \theta_0))>0, m>0より\sin(\theta - \theta_0)=0のとき、つまり\theta = \theta_0, \theta = \theta_0 + \pi\ (\textrm{mod}\ 2\pi)のとき。
\theta = \theta_0では
\frac{\partial^{2} p}{\partial \theta^{2}}=\frac{1}{2 \pi I_{0}(m)} \exp \left(m \cos \left(\theta-\theta_{0}\right)\right)(-m)<0
なので、\theta = \theta_0がP(\theta \mid \theta_0,m)の極大点。
\theta = \theta_0 + \piでは
\frac{\partial^2 p}{\partial \theta^2} = \frac{1}{2\pi I_0(m)}\exp (m \cos (\theta - \theta_0))
m>0より\theta = \theta_0 + \pi\ (\mathrm{mod}\ 2\pi)がP(\theta \mid \theta_0, m)の極小点となる。
演習 2.55
\bar{x}_{1}=\bar{r} \cos \bar{\theta}=\frac{1}{N} \sum_{n=1}^{N} \cos \theta_{n}, \quad \bar{x}_{2}=\bar{r} \sin \bar{\theta}=\frac{1}{N} \sum_{n=1}^{N} \sin \theta_{n} \tag{2.168}
の結果を,
\theta_{0}^{\mathrm{ML}}=\tan ^{-1}\left\{\frac{\sum_{n} \sin \theta_{n}}{\sum_{n} \cos \theta_{n}}\right\} \tag{2.184}
と三角関数の公式
\cos A \cos B+\sin A \sin B=\cos (A-B) \tag{2.178}
と共に用いて,フォン・ミーゼス分布の集中度の最尤推定解m_{\mathrm{ML}}が,A(m_{\mathrm{ML}})=\bar{r}を満たすことを示せ.ただし,\bar{r}は,図2.17のように,2次元ユークリッド平面中の単位ベクトルによって表した観測値の平均の半径である.
(2.187)より
\begin{aligned} A(m_{\mathrm{ML}}) &= \left( \frac{1}{N} \sum_{n=1}^{N} \cos \theta_n \right) \cos \theta_{0}^{\mathrm{ML}} + \left( \frac{1}{N} \sum_{n=1}^{N} \sin \theta_n \right) \sin \theta_{0}^{\mathrm{ML}} \\&= \bar{r} \cos \bar{\theta} \cos \theta_{0}^{\mathrm{ML}} + \bar{r} \sin \bar{\theta} \sin \theta_{0}^{\mathrm{ML}} \\&= \bar{r} \cos \theta_{0}^{\mathrm{ML}} \cos \theta_{0}^{\mathrm{ML}} + \bar{r} \sin \theta_{0}^{\mathrm{ML}} \sin \theta_{0}^{\mathrm{ML}} \\&= \bar{r} \end{aligned}
演習 2.56
ベータ分布
\operatorname{Beta}(\mu \mid a, b)=\frac{\Gamma(a+b)}{\Gamma(a) \Gamma(b)} \mu^{a-1}(1-\mu)^{b-1} \tag{2.13}
ガンマ分布
\operatorname{Gam}(\lambda \mid a, b)=\frac{1}{\Gamma(a)} b^{a} \lambda^{a-1} \exp (-b \lambda) \tag{2.146}
およびフォン・ミーゼス分布
p\left(\theta \mid \theta_{0}, m\right)=\frac{1}{2 \pi I_{0}(m)} \exp \left\{m \cos \left(\theta-\theta_{0}\right)\right\} \tag{2.179}
を指数型分布族の形
p(\mathbf{x} \mid \boldsymbol{\eta})=h(\mathbf{x}) g(\boldsymbol{\eta}) \exp \left\{\boldsymbol{\eta}^{\mathbf{T}} \mathbf{u}(\mathbf{x})\right\} \tag{2.194}
に変形し,これらの分布の自然パラメータを求めよ.
※ 与えられた各分布の形を無理やり指数型分布族の一般形に変形していけば求まる。パラメータの変数名は各分布にしたがって適切に置き換える。
まずベータ分布は
\begin{aligned}
\frac{\Gamma(a+b)}{\Gamma(a) \Gamma(b)} \mu^{a-1}(1-\mu)^{b-1}
&= \frac{\Gamma(a+b)}{\Gamma(a) \Gamma(b)}\exp\left\{ (a-1)\ln \mu +(b-1) \ln(1-\mu) \right\}
\\
&= \frac{\Gamma(a+b)}{\Gamma(a) \Gamma(b)}\exp\left\{ \left(\begin{array}{c}a-1 \\ b-1\end{array}\right)^{\mathrm{T}} \left(\begin{array}{c}\ln \mu \\ \ln (1-\mu)\end{array}\right) \right\}
\end{aligned}
これより、一般形と照らし合わせると、\mathbf{x}\to\mu, \eta \to a,bとして
h(\mu) = 1,\ g(a,b) = \frac{\Gamma(a+b)}{\Gamma(a) \Gamma(b)},\ \boldsymbol{\eta}(a,b) = \left(\begin{array}{c}a-1 \\ b-1\end{array}\right),\ \mathbf{u}(\mu) = \left(\begin{array}{c}\ln \mu \\ \ln (1-\mu)\end{array}\right)
同様にしてガンマ分布は
\begin{aligned}
\frac{1}{\Gamma(a)} b^{a} \lambda^{a-1} \exp (-b \lambda) &= \frac{1}{\Gamma(a)} b^{a} \exp \{(a-1)\ln\lambda - b\lambda\} \\
&= \frac{b^a}{\Gamma(a)} \exp \left\{ \begin{pmatrix}a-1 \\ -b \end{pmatrix}^{\mathrm{T}} \begin{pmatrix}\ln \lambda \\ \lambda \end{pmatrix}\right\}
\end{aligned}
これより
h(\lambda) = 1,\ g(a,b) = \frac{b^a}{\Gamma(a)},\ \boldsymbol{\eta}(a,b) = \left(\begin{array}{c}a-1 \\ -b\end{array}\right),\ \mathbf{u}(\lambda) = \left(\begin{array}{c}\ln \lambda \\ \lambda \end{array}\right)
フォン・ミーゼス分布は
\begin{aligned}
\frac{1}{2 \pi I_{0}(m)} \exp \left\{m \cos \left(\theta-\theta_{0}\right)\right\}
&= \frac{1}{2 \pi I_{0}(m)} \exp \{ m\cos \theta \cos\theta_0 + m\sin\theta \sin\theta_0 \} \\
&= \frac{1}{2 \pi I_{0}(m)} \exp \left\{ \left(\begin{array}{c}m\cos\theta_0 \\ m\sin\theta_0 \end{array}\right)^{\mathrm{T}} \left(\begin{array}{c}\cos\theta \\ \sin\theta \end{array}\right) \right\}
\end{aligned}
これより、
h(\theta) = 1,\ g(\theta_0, m) = \frac{1}{2\pi I_0(m)},\ \boldsymbol{\eta}(\theta_0, m) = \left(\begin{array}{c}m\cos\theta_0 \\ m\sin\theta_0 \end{array}\right),\ \mathbf{u}(\theta) = \left(\begin{array}{c}\cos\theta \\ \sin\theta \end{array}\right)
演習 2.57
多変量ガウス分布は,指数型分布族の形式
p(\mathbf{x} \mid \boldsymbol{\eta})=h(\mathbf{x}) g(\boldsymbol{\eta}) \exp \left\{\boldsymbol{\eta}^{\mathbf{T}} \mathbf{u}(\mathbf{x})\right\} \tag{2.194}
に変形できることを示し,(2.220)–(2.223)と同様に,\boldsymbol{\eta}, \mathbf{u}(\mathbf{x}), h(\mathbf{x})およびg(\boldsymbol{\eta})の式を導出せよ.
多変量ガウス分布の式
\begin{aligned}
N(\mathbf{x} \mid \boldsymbol{\mu}, \mathbf{\Sigma})
&=\frac{1}{(2 \pi)^{D / 2}} \frac{1}{|\mathbf{\Sigma}|^{1 / 2}} \exp \left\{-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu}) \mathbf{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})^{\mathrm{T}}\right\} \\
&=\frac{1}{(2 \pi)^{D / 2}} \frac{1}{|\mathbf{\Sigma}|^{1 / 2}} \exp\left\{-\frac{1}{2} \boldsymbol{\mu}^{\mathrm{T}} \mathbf{\Sigma}^{-1} \boldsymbol{\mu}\right\} \exp \left\{-\frac{1}{2} \mathbf{x}^{\mathrm{T}} \mathbf{\Sigma}^{-1} \mathbf{x}+\mathbf{x}^{\mathrm{T}} \mathbf{\Sigma}^{-1} \boldsymbol{\mu}\right\}
\\
&= \frac{1}{(2 \pi)^{D / 2}} \frac{1}{|\mathbf{\Sigma}|^{1 / 2}} \exp\left\{-\frac{1}{2} \boldsymbol{\mu}^{\mathrm{T}} \mathbf{\Sigma}^{-1} \boldsymbol{\mu}\right\} \exp \left\{-\frac{1}{2}\operatorname{Tr}[\mathbf{xx}^{\mathrm{T}}\mathbf{\Sigma}^{-1}] + \boldsymbol{\mu}^{\mathrm{T}}\mathbf{\Sigma}^{-1}\mathbf{x}\right\}
\end{aligned}
からはじめて、指数型分布族の一般形と比較すると、
h(\mathbf{x}) = (2\pi)^{-D/2}, g(\boldsymbol{\eta}) = |\mathbf{\Sigma}|^{-1 / 2}\exp\left\{-\frac{1}{2} \boldsymbol{\mu}^{\mathrm{T}} \mathbf{\Sigma}^{-1} \boldsymbol{\mu}\right\}
となる。(g(\boldsymbol{\eta})が\boldsymbol{\eta}の関数になっていないことについては後述。)
指数部分について\boldsymbol{\eta}(\boldsymbol{\mu}, \mathbf{\Sigma})と\mathbf{u}(\mathbf{x})の形に分離できるようにすることを考える。
まず\boldsymbol{\mu}^{\mathrm{T}}\mathbf{\Sigma}^{-1}\mathbf{x}を作るため、
\boldsymbol{\eta}(\mu, \mathbf{\Sigma})^{\mathrm{T}} =
\begin{pmatrix}
\mathbf{A} \\
\mathbf{\Sigma}^{-1}\boldsymbol{\mu}
\end{pmatrix}^{\mathrm{T}},\ \mathbf{u}(\mathbf{x}) =
\begin{pmatrix}
\mathbf{B} \\
\mathbf{x}
\end{pmatrix}
とすれば(区分行列(または分割行列、ブロック行列とも)\mathbf{A}, \mathbf{B}の上下の位置は逆でもよい)、\boldsymbol{\eta}^{\mathrm{T}}\mathbf{u} = \mathbf{A}^{\mathrm{T}}\mathbf{B}+\boldsymbol{\mu}^{\mathrm{T}}\mathbf{\Sigma}^{-1}\mathbf{x}となるので、\expの第2項ができる(\because \mathbf{\Sigma}^{-1} = (\mathbf{\Sigma}^{-1})^{\mathrm{T}})。
問題は\mathbf{A}^{\mathrm{T}}\mathbf{B} = \operatorname{Tr}[\mathbf{xx}^{\mathrm{T}}\mathbf{\Sigma}^{-1}]となるような区分行列\mathbf{A}, \mathbf{B}を求めることである。
ここで\mathbf{xx}^{\mathrm{T}} = \mathbf{M}とおくと
\begin{aligned}
\operatorname{Tr}[\mathbf{xx}^{\mathrm{T}}\mathbf{\Sigma}^{-1}] &= \sum_{j=1}^{D}\mathbf{M}_j\mathbf{\Sigma}^{-1}_j \\
&= \begin{pmatrix}\mathbf{m}_{1}^{\mathrm{T}} \mathbf{m}_{2}^{\mathrm{T}} \ldots \mathbf{m}_{n}^{\mathrm{T}}\end{pmatrix} \begin{pmatrix}\boldsymbol{\sigma}_{1}^{-1} \\ \boldsymbol{\sigma}_{2}^{-1} \\ \vdots \\ \boldsymbol{\sigma}_{n}^{-1}\end{pmatrix} \\
&= \begin{pmatrix}\mathbf{m}_{1} \\ \mathbf{m}_{2} \\ \vdots \\ \mathbf{m}_{n}\end{pmatrix}^{\mathrm{T}}\begin{pmatrix}\boldsymbol{\sigma}_{1}^{-1} \\ \boldsymbol{\sigma}_{2}^{-1} \\ \vdots \\ \boldsymbol{\sigma}_{n}^{-1}\end{pmatrix} \\
&= (\operatorname{vec}(\mathbf{M}))^{\mathrm{T}}(\operatorname{vec}(\mathbf{\mathbf{\Sigma}^{-1}})) \\
&= (\operatorname{vec}(\mathbf{xx}^{\mathrm{T}}))^{\mathrm{T}}(\operatorname{vec}(\mathbf{\mathbf{\Sigma}^{-1}}))
\end{aligned}
ここで\operatorname{vec}(\mathbf{A})は任意のm\times n行列\mathbf{A}について要素\mathbf{A} = \{ a_{ij} \}をmn次元列ベクトルに再配列したものである(※vec作用素と呼ばれる。統計のための行列代数下巻 16.2と定理16.2.2を参照。)
たとえば, m = 2, n = 3 の行列\mathbf{A}にvec作用素を適用すると,
\operatorname{vec}(\mathbf{A}) = \begin{pmatrix} \mathbf{a}_{1} \\ \mathbf{a}_{2} \\ \vdots \\ \mathbf{a}_{n} \end{pmatrix}=\begin{pmatrix}a_{11} \\ a_{21} \\ a_{12} \\ a_{22} \\ a_{13} \\ a_{23}\end{pmatrix}
である
これを用いると、求める\boldsymbol{\eta}(\boldsymbol{\mu}, \mathbf{\Sigma}), \mathbf{u}(\mathbf{x})は
\boldsymbol{\eta}(\boldsymbol{\mu}, \mathbf{\Sigma}) =
\begin{pmatrix}
\operatorname{vec}(\mathbf{\Sigma}^{-1}) \\
\mathbf{\Sigma}^{-1}\boldsymbol{\mu}
\end{pmatrix},\ \mathbf{u}(\mathbf{x}) =
\begin{pmatrix}
-\frac{1}{2}\operatorname{vec}(\mathbf{xx}^{\mathrm{T}}) \\
\mathbf{x}
\end{pmatrix}
となる。
(g(\boldsymbol{\eta})が\boldsymbol{\eta}についての式になっていないけれど大丈夫か?という疑問が残るが、vecを作用させる前の\mathbf{\Sigma}^{-1}を\eta_1,\mathbf{\Sigma}^{-1}\boldsymbol{\mu}を\eta_2とすれば, \eta_1, \eta_2の関数としてg(\boldsymbol{\eta})を表現することは可能なのでセーフという理屈らしい……)
演習 2.58
-\nabla \ln g(\boldsymbol{\eta})=\mathbb{E}[\mathbf{u}(\mathbf{x})] \tag{2.226}
は,指数型分布族では\ln g(\boldsymbol{\eta})の負の勾配が,\mathbf{u}(\mathbf{x})の期待値になることを示している.
g(\boldsymbol{\eta}) \int h(\mathbf{x}) \exp \left\{\boldsymbol{\eta}^{\mathrm{T}} \mathbf{u}(\mathbf{x})\right\} \mathrm{d} \mathbf{x}=1 \tag{2.195}
の2階微分を取ることで
-\nabla \nabla \ln g(\boldsymbol{\eta})=\mathbb{E}\left[\mathbf{u}(\mathbf{x}) \mathbf{u}(\mathbf{x})^{\mathrm{T}}\right]-\mathbb{E}[\mathbf{u}(\mathbf{x})] \mathbb{E}\left[\mathbf{u}(\mathbf{x})^{\mathrm{T}}\right]=\operatorname{cov}[\mathbf{u}(\mathbf{x})] \tag{2.300}
を示せ.
※ 先にスカラーについてのベクトルでの微分、ベクトルについてのベクトルでの微分のやり方を知っておく必要がある。
g(\boldsymbol{\eta}) \int h(\mathbf{x}) \exp \left\{\boldsymbol{\eta}^{\mathrm{T}}\mathbf{u}(\mathbf{x})\right\}\mathrm{d}\mathbf{x} = 1なのでg(\boldsymbol{\eta}) \int h(\mathbf{x}) \exp \left\{\boldsymbol{\eta}^{\mathrm{T}}\mathbf{u}(\mathbf{x})\right\}\mathbf{u}(\mathbf{x})はスカラーである。一方で\mathbf{u}(\mathbf{x})は列ベクトルである。\boldsymbol{\eta}でそれぞれを微分する場合、得られる結果はそれぞれ列ベクトルと行列になることに注意する。
スカラーをベクトルで微分する場合は
\frac{\partial f(\mathbf{x})}{\partial \mathbf{x}}=
\begin{pmatrix}
\frac{\partial f}{\partial x_{1}} \\ \vdots \\ \frac{\partial f}{\partial x_{D}}
\end{pmatrix}
これは便宜的に以下のように考えることができる。
\frac{\partial}{\partial \mathbf{x}} f(\mathbf{x})=
\begin{pmatrix}
\frac{\partial}{\partial x_{1}} \\ \vdots \\ \frac{\partial}{\partial x_{D}}
\end{pmatrix} f(\mathbf{x})
これより、
\begin{aligned}
\frac{\partial}{\partial \boldsymbol{\eta}}\underbrace{\exp (\boldsymbol{\eta}^{\mathrm{T}}\mathbf{u}(\mathbf{x}))}_{\textrm{scalar}}
&= \begin{pmatrix}
\frac{\partial}{\partial \eta_{1}} \\ \vdots \\ \frac{\partial}{\partial \eta_{D}}
\end{pmatrix} \exp \left( \sum_{i=1}^{D}\eta_{i}u_i(\mathbf{x}) \right) \\
&= \begin{pmatrix}
\exp \left( \sum_{i=1}^{D}\eta_{i}u_i(\mathbf{x}) \right) u_1(\mathbf{x}) \\ \vdots \\ \exp \left( \sum_{i=1}^{D}\eta_{i}u_i(\mathbf{x}) \right) u_D(\mathbf{x})
\end{pmatrix} \\
&= \exp (\boldsymbol{\eta}^{\mathrm{T}}\mathbf{u}(\mathbf{x})) \mathbf{u}(\mathbf{x})
\end{aligned}
となる。また、ベクトルをベクトルで微分する場合は、教科書によって微分する変数側を行ベクトルとするか、微分される関数側を行ベクトルとするか2通りの表現があるが、ここでは変数側を行ベクトルとする。微分される関数を\mathbf{f}(\mathbf{x})とすると
\frac{\partial}{\partial \mathbf{x}} \mathbf{f}(\mathbf{x})=
\begin{pmatrix}
\frac{\partial f_1}{\partial x_{1}} & \cdots &\frac{\partial f_1}{\partial x_{D}} \\ \vdots & & \vdots \\ \frac{\partial f_D}{\partial x_1} & \ldots & \frac{\partial f_D}{\partial x_D}
\end{pmatrix}
これは便宜的に以下のように考えることができる。
\frac{\partial}{\partial \mathbf{x}} \mathbf{f}(\mathbf{x})=
\begin{pmatrix}
\frac{\partial}{\partial x_{1}} \\ \vdots \\ \frac{\partial}{\partial x_{D}}
\end{pmatrix}
\begin{pmatrix}
f_1(\mathbf{x}) & \ldots & f_D(\mathbf{x})
\end{pmatrix}
これより
\begin{aligned}
\frac{\partial}{\partial \boldsymbol{\eta}}\underbrace{\exp (\boldsymbol{\eta}^{\mathrm{T}}\mathbf{u}(\mathbf{x}))\mathbf{u}(\mathbf{x})}_{\textrm{vector}}
&= \begin{pmatrix}
\frac{\partial}{\partial \eta_{1}} \\ \vdots \\ \frac{\partial}{\partial \eta_{D}}
\end{pmatrix}
\begin{pmatrix}
\exp (\boldsymbol{\eta}^{\mathrm{T}}\mathbf{u}(\mathbf{x})) u_1(\mathbf{x}) & \ldots & \exp (\boldsymbol{\eta}^{\mathrm{T}}\mathbf{u}(\mathbf{x})) u_D(\mathbf{x})
\end{pmatrix} \\
&= \begin{pmatrix}
\exp (\boldsymbol{\eta}^{\mathrm{T}}\mathbf{u}(\mathbf{x})) u_1(\mathbf{x}) u_1(\mathbf{x}) & \ldots & \exp (\boldsymbol{\eta}^{\mathrm{T}}\mathbf{u}(\mathbf{x})) u_D(\mathbf{x}) u_1(\mathbf{x})
\\ \vdots & \ddots & \vdots \\
\exp (\boldsymbol{\eta}^{\mathrm{T}}\mathbf{u}(\mathbf{x})) u_1(\mathbf{x}) u_D(\mathbf{x}) & \ldots & \exp (\boldsymbol{\eta}^{\mathrm{T}}\mathbf{u}(\mathbf{x})) u_D(\mathbf{x}) u_D(\mathbf{x})
\end{pmatrix} \\
&= \exp (\boldsymbol{\eta}^{\mathrm{T}}\mathbf{u}(\mathbf{x})) \mathbf{u}(\mathbf{x}) \mathbf{u}(\mathbf{x})^{\mathrm{T}}
\end{aligned}
となる。
※ PRML P.110の(2.194), (2.195)を用いて(2.224)-(2.226)までの式変形をまず検証する。
(2.195)の両辺を\boldsymbol{\eta}で微分すると
\begin{aligned}
\nabla g(\boldsymbol{\eta}) \int h(\mathbf{x}) \exp \left\{\boldsymbol{\eta}^{\mathrm{T}}\mathbf{u}(\mathbf{x})\right\} \mathrm{d}\mathbf{x}
&+g(\boldsymbol{\eta}) \int h(\mathbf{x}) \exp \left\{\boldsymbol{\eta}^{\mathrm{T}}\mathbf{u}(\mathbf{x})\right\}\mathbf{u}(\mathbf{x}) \mathrm{d}\mathbf{x} = 0 \\
\nabla g(\boldsymbol{\eta}) \int h(\mathbf{x}) \exp \left\{\boldsymbol{\eta}^{\mathrm{T}}\mathbf{u}(\mathbf{x})\right\} \mathrm{d}\mathbf{x} &=
-g(\boldsymbol{\eta}) \int h(\mathbf{x}) \exp \left\{\boldsymbol{\eta}^{\mathrm{T}}\mathbf{u}(\mathbf{x})\right\}\mathbf{u}(\mathbf{x}) \mathrm{d}\mathbf{x} \\
-\frac{\nabla g(\boldsymbol{\eta})}{g(\boldsymbol{\eta})} &= \frac{\int h(\mathbf{x}) \exp \left\{\boldsymbol{\eta}^{\mathrm{T}}\mathbf{u}(\mathbf{x})\right\}\mathbf{u}(\mathbf{x}) \mathrm{d}\mathbf{x}}{\int h(\mathbf{x}) \exp \left\{\boldsymbol{\eta}^{\mathrm{T}}\mathbf{u}(\mathbf{x})\right\}\mathrm{d}\mathbf{x}}
\end{aligned}
ここで再度(2.195)式を使うと
-\frac{\nabla g(\boldsymbol{\eta})}{g(\boldsymbol{\eta})} = g(\boldsymbol{\eta})\int h(\mathbf{x}) \exp \left\{\boldsymbol{\eta}^{\mathrm{T}}\mathbf{u}(\mathbf{x})\right\}\mathbf{u}(\mathbf{x}) \mathrm{d}\mathbf{x} = \mathbb{E}[\mathbf{u}(\mathbf{x})]
これより(2.226)式の-\nabla \ln g(\boldsymbol{\eta}) = \mathbb{E}[\mathbf{u}(\mathbf{x})]が得られる。
続いて2階微分を行うと
\begin{aligned}
-\nabla \nabla \ln g(\boldsymbol{\eta}) &=\nabla g(\boldsymbol{\eta}) \int h(\mathbf{x}) \exp \left\{\boldsymbol{\eta}^{\mathrm{T}} \mathbf{u}(\mathbf{x})\right\} \mathbf{u}(\mathbf{x})^{\mathrm{T}} \mathrm{d}\mathbf{x} \\
&+g(\boldsymbol{\eta}) \int h(\mathbf{x}) \exp \left\{\boldsymbol{\eta}^{\mathrm{T}} \mathbf{u}(\mathbf{x})\right\} \mathbf{u}(\mathbf{x}) \mathbf{u}(\mathbf{x})^{\mathrm{T}} \mathrm{d}\mathbf{x} \\
&=-\mathbb{E}[\mathbf{u}(\mathbf{x})] g(\boldsymbol{\eta}) \int h(\mathbf{x}) \exp \left\{\boldsymbol{\eta}^{\mathrm{T}} \mathbf{u}(\mathbf{x})\right\} \mathbf{u}(\mathbf{x})^{\mathrm{T}} \mathrm{d}\mathbf{x} \\
&+\mathbb{E}\left[\mathbf{u}(\mathbf{x}) \mathbf{u}(\mathbf{x})^{\mathrm{T}}\right] \\ &=\mathbb{E}\left[\mathbf{u}(\mathbf{x}) \mathbf{u}(\mathbf{x})^{\mathrm{T}}\right]-\mathbb{E}[\mathbf{u}(\mathbf{x})]\mathbb{E}\left[\mathbf{u}(\mathbf{x})^{\mathrm{T}}\right] \\
&=\operatorname{cov}[\mathbf{u}(\mathbf{x})]
\end{aligned}
となり、共分散が求まる。
演習 2.59
f(x)が正規化されていれば,密度
p(x|\sigma) = \frac{1}{\sigma}f\left( \frac{x}{\sigma} \right) \tag{2.236}
も正規化されていることを,y=x/\sigmaと変数を変換することで示せ.
\displaystyle \int f(x) dx = 1が成立するときに\displaystyle \int p(x | \sigma) dx = 1であることを示せば良い。
y=x/\sigmaとすると、dy = \frac{1}{\sigma} dxなので
\int p(x|\sigma) dx = \int \sigma \cdot \frac{1}{\sigma}f(y) dy = \int f(y) dy = 1 \hspace{1em} \left( \because \int f(x)dx = 1\right)
したがってp(x|\sigma)が正規化されていることが示された。
演習 2.60
空間\mathbf{x}がいくつかの固定された領域に分割されているとする.このとき,i番目の領域中での密度p(\mathbf{x})は一定の値h_iであり,領域iの体積が\Delta_iであるような,ヒストグラム型の密度モデルを考える.N個の\mathbf{x}の観測値集合があり,領域iに入る観測値が、n_i個であるとする.このとき,密度の正規化制約条件をラグランジュ乗数法によって実現し,\{ h_i \}の最尤推定量の式を導出せよ.
※ ノンパラメトリック法に関連した問題だが、手法として最尤推定を使うので、2.2 多値変数の手法を見直しながらラグランジュ未定乗数法で対数尤度関数の最大化を目指し、\{ h_i \}の最尤推定量を求めていく。しかし、その上で観測データ点\mathbf{x}_nが領域jにあることを示す変数j(n)を使うというのは発想として難しいと思われる。
観測データ点\mathbf{x}_nが領域jにあることを示す変数j(n)を定義すると、点\mathbf{x}_nでの密度p(\mathbf{x})の値はh_{j(n)}で与えられることになる。これを用いると、観測データセット\mathbf{X} = \{ \mathbf{x}_1, \mathbf{x}_2, \ldots, \mathbf{x}_n\}が発生する尤度関数は
\prod_{n=1}^{N}p(\mathbf{x}_n) = \prod_{n=1}^{N}h_{j(n)}
と表せる。すると、この対数尤度関数は
\ln \prod_{n=1}^{N}h_{j(n)} = \sum_{n=1}^{N}\ln h_{j(n)}
となる。
\{h_i\}の最尤推定解を求めるには、制約条件として\displaystyle \sum_{i=1}^{N} h_i \Delta_i = 1となることを考慮しつつ、h_iについての対数尤度関数\displaystyle \sum_{n=1}^{N}\ln h_{j(n)}を最大化する必要がある。すなわちラグランジュ乗数\lambdaを使って
L = \sum_{n=1}^{N}\ln h_{j(n)} + \lambda\left( \sum_{i=1}^{N} h_i \Delta_i - 1 \right)
を最大化する。h_iについての導関数を0とおくと
\frac{\partial L}{\partial h_i} = \frac{\partial}{\partial h_i}\sum_{n=1}^{N}\ln h_{j(n)}+\lambda \Delta_i = 0
ここで\displaystyle \frac{\partial}{\partial h_i}\sum_{n=1}^{N}\ln h_{j(n)}について、j(n)=iとなるような観測値は問題文の設定からn_i個なので
\frac{\partial}{\partial h_i}\sum_{n=1}^{N}\ln h_{j(n)} = n_i \cdot \frac{1}{h_i}
よって\displaystyle \frac{n_i}{h_i}+\lambda \Delta_i = 0。変形するとn_i + \lambda \Delta_i h_i= 0。
iについて和を取ると\sum_i (n_i + \lambda \Delta_i h_i)= N + \lambda = 0となるので、\lambda = -Nが求まる。
以上から、h_iの最尤推定解の式は
h_i = -\frac{n_i}{\lambda \Delta_i} = \frac{n_i}{N\Delta_i}
となる。
演習 2.61
K近傍密度モデルでは,全空間上での積分が発散する変則分布になることを示せ.
【解法1】 ざっくりとしたやり方
あるD次元ユークリッド空間上での未知の確率密度p(\mathbf{x})から観測値の集合が得られていて、この集合からp(\mathbf{x})の値を推定したいとする。(2.246)までの議論から、K近傍法では密度p(\mathbf{x})を推定したい点\mathbf{x}を中心とした半径rのD次元小球の体積V(r)中にK個のデータ点が存在する時(これは言い換えればrはN個の観測値からなるデータセット内で\mathbf{x}から見てK番目に近い点までの距離rをとった時)、密度p(\mathbf{x})の推定値が
p(\mathbf{x}) = \frac{K}{NV(r)}
で与えられることになる。
この問題の題意は「全空間上での積分が発散する」ことを示すことなので、\displaystyle \int p(\mathbf{x})d\mathbf{x} \to \inftyを示せば良い。
ここで、このp(\mathbf{x})を極座標中で考えると、もし十分に大きな半径の値rをとったとき、上式からp(\mathbf{x}) \propto r^{-D}となることがわかる。また、D次元直交座標から極座標に変換するときのヤコビアンは(参考:https://wasan.hatenablog.com/entry/20110321/1300733907)
\begin{aligned}
\left|\frac{\partial\left(x_{1}, \cdots, x_{D}\right)}{\partial\left(r, \theta_{1}, \cdots, \theta_{D-1}\right)}\right| &=\left|\frac{\partial\left(x_{1}, \cdots, x_{D}\right)}{\partial\left(x_{1}, \rho, \theta_{2}, \cdots, \theta_{D-1}\right)}\right|\left|\frac{\partial\left(x_{1}, \rho, \theta_{2}, \cdots, \theta_{D-1}\right)}{\partial\left(r, \theta_{1}, \cdots, \theta_{D-1}\right)}\right| \\ &=\left\{\rho^{D-2} \prod_{i=2}^{D-1}\left(\sin \theta_{i}\right)^{D-i-1}\right\} r \\ &=r\left(r \sin \theta_{1}\right)^{D-2} \prod_{i=2}^{D-1}\left(\sin \theta_{i}\right)^{D-i-1} \quad\left(\because \rho=r \sin \theta_{1}\right) \\ &=r^{D-1} \prod_{i=1}^{D-1}\left(\sin \theta_{i}\right)^{D-i-1} \end{aligned}
より、
d\mathbf{x} = d x_{1} \cdots d x_{D}=r^{D-1}\left\{\prod_{i=1}^{D-1}\left(\sin \theta_{i}\right)^{D-i-1}\right\} d r d \theta_{1} \cdots d \theta_{D-1}
となるので、
\int p(\mathbf{x})d \mathbf{x} \propto \int r^{-D}d\mathbf{x} \propto \int r^{-D}r^{D-1} dr d \theta_{1} \cdots d \theta_{D-1} \propto \int r^{-1}dr = \ln r
これはr\to \inftyの無限大で\ln r \to \infty となるため全空間上での積分は発散することが示された。すなわち、K近傍密度モデルは変則分布であることが示された。
【解法2】より厳密な解き方
https://qiita.com/r-takahama/items/cc03cffd49ad4732ebd3 に書かれてあります。
Discussion
DR YOSHITAKAさん
演習2.60について質問がございます。
ヒストグラムの領域の分割の個数が観測値の個数N個と一致しているのはなぜでしょうか??
M個の分割と考え、制約条件がi=1からM(領域の分割の個数)にするほうが適切に感じます。
最終的な変形にも影響をしますが、結果は同じになると思います。
勉強中の身ですので、間違った質問であれば大変申し訳ございません。