ブログ

当記事は「統計学実践ワークブック(学術図書出版社)」の読解サポートにあたってChapter.$31$の「ベイズ法」に関して演習問題を中心に解説を行います。事前分布を導入し、事後分布や最大事後確率を考えるベイズ的な手法はよく出てくるので、演習を通して抑えておくと良いと思います。

日本統計学会公式認定 統計検定準1級対応 統計学実践ワークブック

3,080円(11/09 13:11時点)

Amazon

本章のまとめ

ベイズ推定量・MAP推定量

具体的な計算を下記で取り扱った。
https://www.hello-statisticians.com/python/stat_program7.html

共役事前分布

「二項分布-ベータ分布」、「正規分布-正規分布」、「ポアソン分布-ガンマ分布」などの共役事前分布などに関して下記で詳しく取り扱った。
https://www.hello-statisticians.com/explain-terms-cat/conjugate_dist1.html

ベイズ線形回帰と予測分布

計算が複雑なのでおそらく準$1$級では出題されないが、応用上重要な結果であるので、結果は抑えておくとよい。
https://www.hello-statisticians.com/explain-terms-cat/bayes_linear_reg1.html

演習問題解説

例31.1

$20$世帯中$x$世帯が保有していたことより、尤度関数$P(x|p)$は下記のように表すことができる。
$$
\large
\begin{align}
P(x|p) = {}_{20} C_{x} p^{x} (1-p)^{20-x}
\end{align}
$$

また、ここで$Be(2,20)$の確率密度関数の$f(p)$は下記のように表される。
$$
\large
\begin{align}
f(p) &= \frac{1}{B(2,20)} p^{2-1} (1-p)^{20-1} \\
&= \frac{1}{B(2,20)} p^{1} (1-p)^{19}
\end{align}
$$

ここまでの議論により、事後分布の確率密度関数$f(p|x)$に関して下記のように計算できる。
$$
\large
\begin{align}
f(p|x) & \propto P(x|p) \times f(p) \\
&= {}_{20} C_{x} p^{x} (1-p)^{20-x} \times \frac{1}{B(2,20)} p^{1} (1-p)^{19} \\
& \propto p^{x+2-1} (1-p)^{40-x-1}
\end{align}
$$

上記より事後分布が$Be(x+2,40-x)$であることがわかる。

例31.2

$1)$
ベータ分布$Be(a,b)$に対して、$\displaystyle E[X]=\frac{a}{a+b}$が成立するので、$Be(2+x,40-x)$に基づくベイズ推定量$\hat{p}_B$は下記のように計算できる。
$$
\large
\begin{align}
\hat{p}_B &= \frac{2+x}{2+x+40-x} \\
&= \frac{2+x}{42}
\end{align}
$$

$2)$
事前分布$Be(2,20)$の期待値は下記のように計算できる。
$$
\large
\begin{align}
\frac{2}{2+20} = \frac{1}{11}
\end{align}
$$

ここでベイズ推定量$\hat{p}_B$は下記のように変形することができる。
$$
\large
\begin{align}
\hat{p}_B &= \frac{2+x}{42} \\
&= \frac{20}{42} \times \frac{x}{20} + \frac{22}{42} \times \frac{1}{11}
\end{align}
$$

上記より、「ベイズ推定量」が「最尤推定値」と「事前分布の期待値」の加重平均で表せることがわかる。

$3)$
ベータ分布$Be(a,b)$のモードは$\displaystyle \frac{a-1}{a+b-2}$のように計算できる。MAP推定量$\hat{p}_{MAP}$は$Be(2+x,40-x)$のモードであるので、下記のように計算できる。
$$
\large
\begin{align}
\hat{p}_{MAP} &= \frac{2+x-1}{2+x+40-x-2} \\
&= \frac{1+x}{40}
\end{align}
$$

・参考
ベータ分布の定義と期待値・分散の計算
ベータ分布の増減表とモード_確率分布の演習

例31.3

$1)$
尤度関数$P(x|\lambda)$は下記のように表すことができる。
$$
\large
\begin{align}
P(x|\lambda) = \prod_{i=1}^{5} \frac{\lambda^{x_i} e^{-\lambda}}{x_i!}
\end{align}
$$

また、ここで$Ga(2,1)$の確率密度関数の$f(\lambda)$は下記のように表される。
$$
\large
\begin{align}
f(\lambda) = \frac{1}{\gamma(2)} \lambda^{2-1} e^{-x}
\end{align}
$$

ここまでの議論により、事後分布の確率密度関数$f(\lambda|x)$に関して下記のように計算できる。
$$
\large
\begin{align}
f(\lambda|x) & \propto P(x|\lambda) \times f(\lambda) \\
&= \prod_{i=1}^{5} \frac{\lambda^{x_i} e^{-\lambda}}{x_i!} \times \frac{1}{\gamma(2)} \lambda^{2-1} e^{-x} \\
& \propto \lambda^{\sum_{i=1}^{5} x_i + 2 -1} e^{-x-5x} \\
&= \lambda^{3+7+1+4+5+2-1} e^{-6x} \\
&= \lambda^{22-1} e^{-6x}
\end{align}
$$

上記より事後分布が$Ga(22,1/6)$であることがわかる。

$2)$
ガンマ分布$Be(\nu,\alpha)$のモードは$\displaystyle E[X] = \nu \alpha$のように計算できる。よって、$Ga(22,1/6)$の期待値は下記のように計算できる。
$$
\large
\begin{align}
E[X] &= \nu \alpha \\
&= 22 \times \frac{1}{6} \\
&= \frac{11}{3}
\end{align}
$$

$3)$
ガンマ分布$Ga(22,1/6)$のモードは$\displaystyle (\nu-1) \alpha$のように計算できる。MAP推定量$\hat{\lambda}_{MAP}$は$Be(22,1/6)$のモードであるので、下記のように計算できる。
$$
\large
\begin{align}
\hat{\lambda}_{MAP} &= (22-1) \times \frac{1}{6} \\
&= \frac{7}{2}
\end{align}
$$

・参考
ガンマ分布の定義と期待値・分散の計算
ガンマ分布の増減表とモード_確率分布の演習

問31.1

$1)$
一様分布$U(0,1)$にはベータ分布$Be(1,1)$が対応する。ここで、尤度関数が$L(\theta)={}_{10} C_{3} \theta^{3} (1-\theta)^{7}$なので、パラメータの事後分布を$f(\theta|x)$のようにおくと、$f(\theta|x)$は下記のように計算できる。
$$
\large
\begin{align}
f(\theta|x) & \propto \theta^{3} (1-\theta)^{7} \\
&= \theta^{4-1} (1-\theta)^{8-1}
\end{align}
$$

上記より、事後分布は$Be(4,8)$であることがわかる。また、これに対応するベイズ推定量は$Be(a,b)$の期待値が$\displaystyle E[X]=\frac{a}{a+b}$であることより、$\displaystyle \frac{4}{4+8}=\frac{1}{3}$であることがわかる。

$2)$
$1)$と同様に事後分布は下記のように計算できる。
$$
\large
\begin{align}
f(\theta|x) & \propto \theta^{1-1} (1-\theta)^{4-1} \times \theta^{3} (1-\theta)^{7} \\
&= \theta^{4-1} (1-\theta)^{11-1}
\end{align}
$$

上記位より事後分布は$Be(4,11)$であることがわかる。また、これに対応するMAP推定量は$\displaystyle \frac{4-1}{4+11-2}=\frac{3}{13}$であることがわかる。

・参考
ベータ分布の定義と期待値・分散の計算
ベータ分布の増減表とモード_確率分布の演習

問31.2

$1)$
$\lambda$の事後分布$\pi(\lambda|x)$は下記のように計算できる。
$$
\large
\begin{align}
\pi(\lambda|x) & \propto \prod_{i=1}^{5} \frac{\lambda^{x_i} e^{-\lambda}}{x_i!} \times \lambda^{2-1} e^{-3 \lambda} \\
& \propto \lambda^{5+3+4+1+2+2-1} e^{-3 \lambda – 5 \lambda} \\
&= \lambda^{17-1} e^{-8 \lambda}
\end{align}
$$

よって、$\lambda$の事後分布は$Ga(17,1/8)$であり、確率密度関数$\pi(\lambda|x)$は下記のように表せる。
$$
\large
\begin{align}
\pi(\lambda|x) = \frac{8^{17}}{\gamma(17)} \lambda^{17-1} e^{-8 \lambda}
\end{align}
$$

$2)$
$$
\large
\begin{align}
f(x_6|\lambda) = \frac{\lambda^{x_6} e^{-\lambda}}{x_6!}
\end{align}
$$

上記のように$f(x_6|\lambda)$が表せることからベイズ予測分布$f(x_6|x)$の確率密度関数は下記のように計算できる。
$$
\large
\begin{align}
f(x_6|x) &= \int_{0}^{\infty} f(x_6|\lambda) \pi(\lambda|x) d \lambda \\
&= \int_{0}^{\infty} \frac{\lambda^{x_6} e^{-\lambda}}{x_6!} \times \frac{8^{17}}{\gamma(17)} \lambda^{17-1} e^{-8 \lambda}d \lambda \\
&= \frac{8^{17}}{x_6! \gamma(17)} \int_{0}^{\infty} \lambda^{17+x_6-1} e^{-9\lambda} d \lambda \\
&= \frac{8^{17}}{x_6! \gamma(17)} \times \frac{\gamma(17+x_6)}{9^{17+x_6}} \\
&= \frac{\gamma(17+x_6)}{x_6! \gamma(17)} \cdot \left(\frac{8}{9}\right)^{17} \cdot \left(\frac{1}{9}\right)^{x_6}
\end{align}
$$

上記の計算にあたっては、ガンマ分布$Ga(\nu,\alpha)$の確率密度関数の全区間での積分に関する下記の関係式を用いた。
$$
\large
\begin{align}
\int_{0}^{\infty} \frac{1}{\gamma(\nu) \alpha^{\nu}} x^{\nu-1} e^{-\frac{x}{\alpha}} dx &= 1 \\
\int_{0}^{\infty} x^{\nu-1} e^{-\frac{x}{\alpha}} &= \gamma(\nu) \alpha^{\nu}
\end{align}
$$

問31.3

$1)$
$\alpha(x^{(t)},y)$は下記のように表現できる。
$$
\large
\begin{align}
\alpha(x^{(t)},y) &= \min \left\{ 1, \frac{\pi(y)}{\pi(x^{(t)})} \right\} \\
&= \min \left\{ 1, \frac{\frac{1}{8} \varphi(y+5) + \frac{3}{4} \varphi(y) + \frac{1}{8} \varphi(y-5)}{\frac{1}{8} \varphi(x^{(t)}+5) + \frac{3}{4} \varphi(x^{(t)}) + \frac{1}{8} \varphi(x^{(t)}-5)} \right\}
\end{align}
$$

$2)$
下記のように対応する。
ア) $a=1$
イ) $a=10$
ウ) $a=100$
エ) $a=0.1$

基本的にはステップ幅が大きい方が値の変動が起こりやすいが、$a=100$のように移動が大きい場合は採択確率が下がる場合があることに注意すると良い。

$3)$
繰り返し回数が少ない場合は初期値に左右されることがある。

問31.4

$$
\large
\begin{align}
f(x,y) &= \frac{1}{2 \pi \sqrt{1-0.5^2}} \exp \left( -\frac{x^2 – 2 \cdot 0.5 x(y-2) + (y-2)^2)}{2(1-0.5^2)} \right) \\
f(x) &= \frac{1}{\sqrt{2 \pi}} \exp \left( -\frac{x^2}{2} \right) \\
f(y) &= \frac{1}{\sqrt{2 \pi}} \exp \left( -\frac{(y-2)^2}{2} \right)
\end{align}
$$

また、条件付き確率密度関数$f(x|y), f(y|x)$に関して下記が成立する。
$$
\large
\begin{align}
f(x|y) &= \frac{f(x,y)}{f(y)} \\
&= \frac{1}{\sqrt{2 \pi (1-0.5^2)}} \exp \left( -\frac{x^2 – 2 \cdot 0.5 x(y-2) + (y-2)^2}{2(1-0.5^2)} – \frac{(y-2)^2}{2} \right) \quad (1) \\
f(y|x) &= \frac{f(x,y)}{f(x)} \\
&= \frac{1}{\sqrt{2 \pi (1-0.5^2)}} \exp \left( -\frac{x^2 – 2 \cdot 0.5 x(y-2) + (y-2)^2}{2(1-0.5^2)} – \frac{x^2}{2} \right) \quad (2)
\end{align}
$$

ここまでの式を元に各設問に関して考える。

$1)$
$f(x_{(t+1)}|y_{(t)})$より、$(1)$式を元に考える。$(1)$式は正規分布を表すことより、平均$\mu_{x_{(t+1)}|y_{(t)}}$と分散$\sigma_{x_{(t+1)}|y_{(t)}}^2$をそれぞれ考えれば良い。式の形より、$\sigma_{x_{(t+1)}|y_{(t)}}^2$は下記のように得られる。
$$
\large
\begin{align}
\sigma_{x_{(t+1)}|y_{(t)}}^2 &= 1-0.5^2 \\
&= \frac{3}{4}
\end{align}
$$

次に$\mu_{x_{(t+1)}|y_{(t)}}$を求めるにあたって、$(1)$式の変数を置き換えた式の指数関数の中身を$x_{(t+1)}$に関して整理することを考える。
$$
\large
\begin{align}
& -\frac{x_{(t+1)}^2 – 2 \cdot 0.5 x_{(t+1)}(y_{(t)}-2) + (y_{(t)}-2)^2}{2(1-0.5^2)} – \frac{(y_{(t)}-2)^2}{2} \\
&= -\frac{\left( x_{(t+1)} – \frac{1}{2}(y_{(t)}-2) \right)^2 + …}{2(1-0.5^2)}
\end{align}
$$

上記より、$\mu_{x_{(t+1)}|y_{(t)}}$は下記のように得られる
$$
\large
\begin{align}
\mu_{x_{(t+1)}|y_{(t)}} = \frac{1}{2}(y_{(t)}-2)
\end{align}
$$

$2)$
$f(y_{(t+1)}|x_{(t+1)})$より、$(2)$式を元に考える。$(2)$式は正規分布を表すことより、平均$\mu_{y_{(t+1)}|x_{(t+1)}}$と分散$\sigma_{y_{(t+1)}|x_{(t+1)}}^2$をそれぞれ考えれば良い。式の形より、$\sigma_{y_{(t+1)}|x_{(t+1)}}^2$は下記のように得られる。
$$
\large
\begin{align}
\sigma_{y_{(t+1)}|x_{(t+1)}}^2 &= 1-0.5^2 \\
&= \frac{3}{4}
\end{align}
$$

次に$\mu_{y_{(t+1)}|x_{(t+1)}}$を求めるにあたって、$(2)$式の変数を置き換えた式の指数関数の中身を$x_{(t+1)}$に関して整理することを考える。
$$
\large
\begin{align}
& -\frac{x_{(t+1)}^2 – 2 \cdot 0.5 x_{(t+1)}(y_{(t)}-2) + (y_{(t)}-2)^2}{2(1-0.5^2)} – \frac{x_{(t+1)}^2}{2} \\
& = -\frac{(y_{(t)} – 2 – \frac{1}{2}x_{(t+1)})^2 + …}{2(1-0.5^2)}
\end{align}
$$

上記より、$\mu_{y_{(t+1)}|x_{(t+1)}}$は下記のように得られる
$$
\large
\begin{align}
\mu_{y_{(t+1)}|x_{(t+1)}} = \frac{1}{2}x_{(t+1)} + 2
\end{align}
$$

参考

・準1級関連まとめ
https://www.hello-statisticians.com/toukeikentei-semi1

・Pythonで学ぶ統計モデリング（運営者作成）
https://www.amazon.co.jp/dp/B08FYMTYBW/

・共役事前分布
https://www.hello-statisticians.com/explain-terms-cat/conjugate_dist1.html