当記事は「基礎統計学Ⅲ 自然科学の統計学(東京大学出版会)」の読解サポートにあたってChapter.$4$の「最尤法」の章末問題の解説について行います。
基本的には書籍の購入者向けの解答例・解説なので、まだ入手されていない方は下記より入手をご検討ください。また、解説はあくまでサイト運営者が独自に作成したものであり、書籍の公式ページではないことにご注意ください。
章末の演習問題について
問題4.1の解答例
下記を実行すればよい。
import numpy as np
import matplotlib.pyplot as plt
x = np.array([3.03, 3.09, 3.15, 3.20, 3.25, 3.30, 3.34, 3.37, 3.41, 3.45, 3.48])
p = np.array([0., 0., 0.11, 0.28, 0.48, 0.71, 0.83, 0.91, 0.95, 0.98, 1.])
plt.plot(np.log(x),p)
plt.show()・実行結果
問題4.2の解答例
i)
$i$回目の実験における表が出るまでの回数を$x_i$とする。この確率を$P(X=x_i|p)$とすると、$P(X=x_i|p)=f(x_i)$より下記のように表せる。
$$
\begin{align}
P(X=x_i|p) &= f(x_i) \\
&= p(1-p)^{x_i}
\end{align}
$$
このとき、$n$回の実験に関する同時確率(joint probability)は$P(X=x_1,X=x_2,…,X=x_n|p)$と表すことができ、$X_1$〜$X_n$が独立であるので下記のように式変形ができる。
$$
\begin{align}
P(X=x_1,X=x_2,…,X=x_n|p) &= P(X=x_1|p)P(X=x_2|p)…P(X=x_n|p) \\
&= f(x_1)f(x_2)…f(x_n) \\
&= p^n(1-p)^{\sum_{i=1}^{n} x_i-1} \\
&= p^n(1-p)^{T-n}
\end{align}
$$
上記では$\displaystyle T = \sum_{i=1}^{n} x_i$を利用した。また、このときの尤度を$L(p)$とすると、$L(p)=P(X=x_1,X=x_2,…,X=x_n|p)$なので、$L(p)=p^n(1-p)^{T-n}$となる。以下、対数尤度の微分$\displaystyle \frac{\delta \log{L(p)}}{\delta p} = 0$を満たす$p$を求める。
$$
\begin{align}
\log{L(p)} &= \log{p^n(1-p)^{T-n}} \\
&= n\log{p}+(T-n)\log{(1-p)}
\end{align}
$$
$$
\begin{align}
\frac{\delta \log{L(p)}}{\delta p} &= \frac{n}{p}-\frac{T-n}{1-p} = 0 \\
\frac{n}{p} &= \frac{T-n}{1-p} \\
n(1-p) &= (T-n)p \\
n-np &= Tp-np \\
n &= Tp \\
p &= \frac{n}{T}
\end{align}
$$
上記より、$\displaystyle \tilde{p} = \frac{n}{T}$を示すことができる。
また、ここでは$n=100$、$T=192$より、$\displaystyle \tilde{p} = \frac{100}{192} = 0.52083…$となる。
問題4.3の解答例
i)
平均$\mu=0$、分散$\sigma^2$の正規分布は下記のように書ける。
$$
\begin{align}
P(x|\mu=0,\sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left( -\frac{x^2}{2\sigma^2} \right)
\end{align}
$$
このとき、$P(x_1,x_2,…,x_n|0,\sigma^2)$の同時確率は下記のように計算できる。
$$
\begin{align}
P(x_1,x_2,…,x_n|0,\sigma^2) &= P(x_1|0,\sigma^2)P(x_2|0,\sigma^2)…P(x_n|0,\sigma^2) \\
&= \prod_{i=1}^{n} \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left( -\frac{x_i^2}{2\sigma^2} \right)
\end{align}
$$
上記の同時確率は$\sigma^2$に着目すると尤度と見なせるので、$L(\sigma^2)=P(x_1,x_2,…,x_n|0,\sigma^2)$が成立する。このとき積の形を和の形に直すにあたって、$L(\sigma^2)$の対数の$\log{L(\sigma^2)}$を下記で計算する。
$$
\begin{align}
\log{L(\sigma^2)} &= \log\left( \prod_{i=1}^{n} \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left( -\frac{x_i^2}{2\sigma^2} \right) \right) \\
&= \log\left(\prod_{i=1}^{n} \frac{1}{\sqrt{2\pi\sigma^2}} \right) + \log\left(\prod_{i=1}^{n} \exp\left( -\frac{x_i^2}{2\sigma^2} \right) \right) \\
&= \log\left( \frac{1}{\sqrt{2\pi\sigma^2}} \right)^n + \log\left(\exp\left( \sum_{i=1}^{n} -\frac{x_i^2}{2\sigma^2} \right)\right) \\
&= \log\left( 2\pi\sigma^2 \right)^{-n/2} – \sum_{i=1}^{n} \frac{x_i^2}{2\sigma^2} \\
&= – \frac{n}{2}\log\left( 2\pi \right) – \frac{n}{2}\log\left( \sigma^2 \right) – \sum_{i=1}^{n} \frac{x_i^2}{2\sigma^2}
\end{align}
$$
次に上記で求めた$\log{L(\sigma^2)}$を$\sigma^2$で微分する。
$$
\begin{align}
\frac{\delta \log{L(\sigma^2)}}{\delta \sigma^2} = – \frac{n}{2\sigma^2} + \sum_{i=1}^{n} \frac{x_i^2}{2(\sigma^2)^2}
\end{align}
$$
ここで$\displaystyle \frac{\delta \log{L(\sigma^2)}}{\delta \sigma^2}=0$を考える。
$$
\begin{align}
& \frac{\delta \log{L(\sigma^2)}}{\delta \sigma^2} = 0 \\
& – \frac{n}{2\sigma^2} + \sum_{i=1}^{n} \frac{x_i^2}{2(\sigma^2)^2} = 0 \\
& \sum_{i=1}^{n} \frac{x_i^2}{2(\sigma^2)^2} = \frac{n}{2\sigma^2} \\
& \sum_{i=1}^{n} x_i^2 = \frac{n \times 2(\sigma^2)^2}{2\sigma^2} \\
& \sum_{i=1}^{n} x_i^2 = n\sigma^2 \\
& \sigma^2 = \frac{\sum_{i=1}^{n} x_i^2}{n}
\end{align}
$$
上記より、$\displaystyle \tilde{\sigma} = \sqrt{\frac{\sum_{i=1}^{n} x_i^2}{n}}$を求めることができる。
問題4.4の解答例
i)
指数分布の確率密度関数を$f(t)=P(t|\lambda)$としたとき、下記のように表すことができる。
$$
\begin{align}
P(t|\lambda) = \lambda e^{-\lambda t} \\
(t \geq 0)
\end{align}
$$
この際に観測されたサンプルを$t_1$、$t_2$、…、$t_{10}$、尤度を$L(\lambda)$とすると、尤度はサンプルが観測された同時確率をパラメータの関数と読み替えたものなので下記が成立する。
$$
\begin{align}
L(\lambda) &= P(t_1|\lambda)P(t_2|\lambda)…P(t_{10}|\lambda) \\
&= \prod_{i=1}^{10} \lambda e^{-\lambda t_i}
\end{align}
$$
ここで上記の対数をとった対数尤度の$\log{L(\lambda)}$について計算すると下記のようになる。
$$
\begin{align}
\log{L(\lambda)} &= \log\left( \prod_{i=1}^{10} \lambda e^{-\lambda t_i} \right) \\
&= \log\left( \lambda^{10} \prod_{i=1}^{10}e^{-\lambda t_i} \right) \\
&= \log\left( \lambda^{10} e^{\sum_{i=1}^{10} -\lambda t_i} \right) \\
&= 10\log{\lambda} – \sum_{i=1}^{10} \lambda t_i \\
&= 10\log{\lambda} – \lambda \sum_{i=1}^{10} t_i
\end{align}
$$
上記を$\lambda$について微分し、$0$になる$\lambda$が$\tilde{\lambda}$となる。
$$
\begin{align}
& \frac{\delta \log{L(\lambda)}}{\delta \lambda} = 0 \\
& \frac{10}{\lambda} – \sum_{i=1}^{10} t_i = 0 \\
& \frac{10}{\lambda} = \sum_{i=1}^{10} t_i \\
& \lambda = \frac{10}{\sum_{i=1}^{10} t_i}
\end{align}
$$
上記より、$\displaystyle \tilde{\lambda} = \frac{10}{\sum_{i=1}^{10} t_i}$が成立する。これを元に下記のように$\tilde{\lambda}$が推定できる。
$$
\begin{align}
\tilde{\lambda} &= \frac{10}{523+971+878+130+696+201+1643+491+665+431} \\
&= \frac{10}{6629} \\
&= 0.0015085…
\end{align}
$$
ⅱ)
指数分布の累積分布関数を$F(t)=F(X \leq t)$とすると、下記のように表すことができる。
$$
\begin{align}
F(X \leq t) &= f(t) \\
&= 1-e^{-\lambda t} \\
F(X \geq t) &= 1-F(X \leq t) \\
&= 1-(1-e^{-\lambda t}) \\
&= e^{-\lambda t}
\end{align}
$$
よって、$F(X \geq 1000)$は下記のように表すことができる。
$$
\begin{align}
F(X \geq 1000) &= e^{-(10 \times 1000)/6629} \\
&= e^{-(10000)/6629} \\
&= 0.221236…
\end{align}
$$
問題4.5の解答例
$\tilde{\theta}=\max{x_i}$が最尤推定量となる。
まとめ
Chapter.$4$で取り扱った最尤法は得られたサンプルの背後に仮定した確率分布のパラメータの推定を行う手法です。同時確率を尤度と読み替えるところがなかなかイメージがつかないかもしれませんが、最重要項目なのでこの辺の理解を中心に理解を進めていくと良いと思います。