2015-10-05

2014年度統計検定1級大問2(統計数理)((2)まで)

大問2(問題文)

連続型確率変数 $X$ の確率密度関数が、 $m>0$ として
$\displaystyle f(x)= \!\begin{cases} \frac{1}{\Gamma(m)}x^{m-1}e^{-x} & (x>0) \\ 0 & (x \le 0) \end{cases}$
で与えられるとき、 $X$ はパラメータ $m$ のガンマ分布にしたがうという。ここで $\Gamma(m)$ はガンマ関数を表し、 $m$ が正の整数ならば $\Gamma(m)=(m-1)!$ である。

(1)パラメータ $m$ のガンマ分布のモーメント母関数は、 $t<1$ に対して $M_x(t)=(1-t)^{-m}$ であることを示せ。

(2)パラメータ $m$ のガンマ分布の平均と分散をそれぞれ求めよ。

(3) $n+1$ 個の確率変数 $X_1,...,X_{n+1}$ が、それぞれ互いに独立に $m=1$ のガンマ分布に従うとき、
$T=X_1+...+X_n+X_{n+1}$
および、 $Y_1=X_1/T, Y_2=(X_!+X_2)/T,...,Y_n=(X_1+...+X_n)/T$ とする。このとき、 $T$ の確率分布は何か。
また、(Y_1,...,Y_n)と $T$ とは互いに独立であることを示せ。

(4)上問(3)で定義された $(Y_1,...,Y_n)$ の同時分布は、互いに独立に区間(0,1)上の一様分布に従う確率変数 $U_1,...,U_n$ の順序統計量[tex: U_{(1)}<...

解説

(1)

モーメント母関数の定義 $M_X(t)=E(e^{tX})$ を利用して計算すればよいです。
このモーメント母関数の定義がなぜこうなるのか、というお話ですが、テイラー展開した時のn次の項がn次のモーメントに対応しているようです(説明になっているようではぐらかしている)。なので $t$ には0近傍で求めるという以外にあまり意味はないです。詳しくはこの辺(http://www.aandt.co.jp/jpn/qc/basic/bokansu.htm)などを参考にして下さい。
説明はこの辺で終わりにして愚直に計算します。
${\displaystyle \begin{eqnarray*} M_X(t) &=&\int_0^\infty \frac{1}{\Gamma(m)}x^{m-1}e^{-x} e^{tx}dx\\ &=&\int_0^\infty \frac{1}{\Gamma(m)}x^{m-1}e^{-(1-t)x}dx\\ \end{eqnarray*} }$
ここで $(1-t)x=y$ と置換すると
${\displaystyle \begin{eqnarray*} \int_0^\infty \frac{1}{\Gamma(m)}x^{m-1}e^{-(1-t)x}dx &=&\int_0^\infty \frac{1}{\Gamma(m)} \left(\frac{y}{1-t}\right)^{m-1} e^{-y}\frac{dy}{1-t}\\ &=&\frac{1}{(1-y)^m}\color{red}{\int_0^\infty \frac{1}{\Gamma(m)} y^{m-1} e^{-y}dy}\\ &=&(1-y)^{-m} \end{eqnarray*} }$
赤い部分はガンマ分布そのものを表しますから、全区間での積分結果は1になります。
したがってパラメータ $m$ のガンマ分布のモーメント母関数は、 $t<1$ に対して $M_x(t)=(1-t)^{-m}$ と示せます。

(2)

モーメント母関数をtについて微分して求めるのが一番楽だと思いますが、別に平均と分散の定義通りに求めても良いと思います。両方書いときますね。

(i)モーメント母関数を利用する場合

平均 $E(X)$ はモーメント母関数をtについて微分して $t=0$ における傾きを出せばよく(1次のモーメント)、
$E(X)=M_X'(0)=(-m)*(-1)*(1-0)^{-m-1}=m$
分散 $V(X)$ は分散の定義より $V(X)=E(X^2)-E(X)^2$ から、 $E(X^2)$ を求めればよく、これはモーメント母関数のtに関する2階微分(2次のモーメント)なので
$E(X^2)=M_X''(0)=(-m)*(-m-1)*(-1)*(-1)*(1-0)^{-m-2}=m(m+1)$
したがって、 $V(X)=E(X^2)-E(X)^2=m(m+1)-m^2=m$
よって平均・分散ともに $m$ と求まります。

(ii)平均と分散の定義から直接求める場合

平均の定義より $E(X)=\inf_0^\infty x f(x)dx (f(x)は確率密度関数)$
これにガンマ分布を代入して
${\displaystyle \begin{eqnarray*} E(X) &=&\inf_0^\infty x f(x)dx\\ &=&\inf_0^\infty x \frac{1}{\Gamma(m)}x^{m-1}e^{-x} dx\\ &=&\inf_0^\infty \frac{1}{\Gamma(m)}x^{m}e^{-x} dx\\ &=&\inf_0^\infty \frac{m}{\Gamma(m+1)}x^{m}e^{-x} dx\\ &=&m \color{red}{\inf_0^\infty \frac{1}{\Gamma(m+1)}x^{m}e^{-x} dx}\\ &=&m \end{eqnarray*} }$
赤字の部分はパラメータm+1のガンマ分布の非零である全区間における積分なので1になり、このように平均 $E(X)$ を求めることができます。
分散は $E(X^2)$ を求めたいので
${\displaystyle \begin{eqnarray*} E(X) &=&\inf_0^\infty x^2 f(x)dx\\ &=&\inf_0^\infty x^2 \frac{1}{\Gamma(m)}x^{m-1}e^{-x} dx\\ &=&\inf_0^\infty \frac{1}{\Gamma(m)}x^{m+1}e^{-x} dx\\ &=&\inf_0^\infty \frac{m(m+1)}{\Gamma(m+2)}x^{m+1}e^{-x} dx\\ &=&m(m+1) \color{red}{\inf_0^\infty \frac{1}{\Gamma(m+1)}x^{m}e^{-x} dx}\\ &=&m(m+1) \end{eqnarray*} }$
赤字の部分はパラメータm+2のガンマ分布の非零である全区間における積分なので1になり、このように $E(X^2)$ を求めることができます。
あとは $V(X)=E(X^2)-E(X)^2=m(m+1)-m^2=m$ と求まります。

(3)以降は明日時間を見つけてやります。

2015-10-04

2014年度統計検定1級大問1(統計数理)

統計検定

注)僕が今年受けるかどうかは未定です。

ぼちぼちkaggleをやりながら統計検定1級の勉強をしているのですが、日本統計学会の出している昨年度の略解(
http://www.toukei-kentei.jp/about/pastpaper/2014n/ans2014n_grade1_suri.pdf
)のあまりの不親切さと未だに刊行されない2014年度版過去問集に憤りを覚えたのでまとめてみます。
なお、以下のブログを大いに参考にしています。blog.goo.ne.jp

先駆者がいらっしゃることですし、本ブログの方針に則り、できるだけ詳しく展開を書くことで差別化を図ろうと思います。

おことわり
・勉強がてら調べつつ解いていくつもりなので内容の正確さはかなり微妙ですのでご了承下さい。間違いに対するご指摘・お叱りは歓迎致します。
・2014年度分が解説付きで過去問集に掲載された後に関してですが、著作権等で問題があるようならばこれらの記事を削除します。(試験問題とその解説ってどういう扱いなんでしょうかね？これがだめならPRMLの方もダメな気がする。。。)

大問1(問題)

区間(0,1)上の一様分布に従う互いに独立な確率変数 $U, V, W$ に対し、 $\alpha, \beta, \gamma$ を正の定数として、
$X=U^\alpha, Y=V^\beta, Z=W^\gamma$ とする。このとき以下の設問に答えよ

(1) $\alpha=2,\beta=3$ のとき、 $U=u$ が与えられた下での条件付き確率 $P(X>Y|U=u)$ を求めよ。
これを用いて $X>Y$ となる確率 $P(X>Y)$ を求めよ。

(2)一般の $\alpha, \beta, \gamma$ に対し、 $X$ が $X, Y, Z$ の中で最大となる確率 $P(X=max(X,Y,Z))$ を求めよ。

(3)上問(2)を拡張し、 $U_i(i = 1,...,n)$ を互いに独立に区間(0,1) 上の一様分布に従う $n$ 個の確率変数とし、 $\alpha_i(i = 1,...,n)$ を正の定数として、 $X_i=U_i^{\alpha_i}(i = 1,...,n)$ とするとき、 $X_1$ が $X_1,...,X_n$ の中で最大となる確率 $P(X_1=max(X_1,...,X_n))$ を $\alpha_1,...,\alpha_n$ の関数として求めよ。

解説

(1)

$U=u$ が与えられた下でのVの確率を $v$ とします。
$\alpha=2,\beta=3$ が与えられていますから、 $X>Y$ が成り立つ時、
$u,v$ には $1>u^2>v^3>0$ の関係があります。
ここで、 $u,v$ は確率なので0から1の範囲をとりますから三乗根をとっても不等号は変わらず
$u^{2/3}>v>0$
が成り立ちます。
よって $u$ が与えられた時、vのとりうる範囲は $u^{2/3}>v>0$ [1]。
Vは区間(0,1)上の一様分布なので、vが[1]の範囲にある確率は $\frac{u^{2/3}-0}{1-0}=u^{2/3}$ 。
(もっと詳しく言えば、「始点0、終点1の長さ1の線分上における、0から $u^{2/3}$ までの部分が占める割合」というイメージです)
以上より条件付き確率 $P(X>Y|U=u)$ ( $U=u$ という条件が与えられた下で、確率変数Xが確率変数Yよりも大きくなる確率)は $u^{2/3}$ となります。

あとは $P(X>Y)$ を求めますが、ここまでくれば簡単で、あり得る全ての $u$ の場合の条件付き確率 $P(X>Y|U=u)$ を足し合わせればよいです。つまり積分すればいいですね。 $u$ が確率だったことを思い出すと、積分区間は0から1までで良いとわかります。
$P(X>Y)=\int_0^1P(X>Y|U=u)du=\int_0^1 u^{2/3}du=[\frac{3}{5}u^{5/3} ]_0^1=\frac{3}{5}$
より、 $P(X>Y)=\frac{3}{5}$ とわかります。

(2)

$P(X=max(X,Y,Z))$ となる確率を求める問題ですが、このままでは少なくとも僕には解けないので $P(X>max(Y,Z))$ と読み替えます。
(ここの部分は2014年度、統計数理、問１の解説 - the BLOG for 統計検定1級を読んでようやく理解しました)
(左の式は「X,Y,Zの中でXが最大となる確率」、右の式は「XがYとZの最大値よりも大きい確率」という意味です。どちらも同じですよね？)

(1)に倣って、 $U=u$ が与えられた下での条件付き確率 $P(X>max(Y,Z)|U=u)$ を求めて積分する方針でいきましょう。
$U=u$ が与えられた時の $V,W$ の確率をそれぞれ $v,w$ とすると、
$u^\alpha>max(v^\beta,w^\gamma)$ となる $v,w$ の条件を調べれば良さそうです。
maxが邪魔なので、ひとまずこれを取り払って $u^\alpha>v^\beta$ と $u^\alpha>w^\gamma$ に分解して考えましょう。
$u^\alpha>v^\beta$ を満たす確率(すなわち $P(X>Y|U=u)$ )は(1)と同様に、 $u^{\alpha/\beta}>v>0$ から、
$u^{\alpha/\beta}$ [2]。
$u^\alpha>w^\gamma$ を満たす確率(すなわち $P(X>Z|U=u)$ )は(1)と同様に、 $u^{\alpha/\gamma}>w>0$ から、
$u^{\alpha/\gamma}$ [3]。

ここで少し考えると、「 $u^\alpha$ が、 $v^\beta$ と $w^\gamma$ の最大値よりも大きい」ということは「 $u^\alpha$ が、 $u^\alpha>v^\beta$ かつ $u^\alpha>w^\gamma$ を満たす」ということと同じであると気づきます。
したがって、 $u^\alpha>max(v^\beta,w^\gamma)$ となる確率(すなわち $P(X>max(Y,Z)|U=u)$ )は、[2],[3]の同時確率より、
$u^{\alpha/\beta}*u^{\alpha/\gamma}=u^{\frac{\alpha\gamma+\alpha\beta}{\beta\gamma}}$
あとは積分して、
${\displaystyle \begin{eqnarray*} P(X>max(Y,Z)) &=&\int_0^1 P(X>max(Y,Z)|U=u)du\\ &=& \int_0^1 u^{\frac{\alpha\gamma+\alpha\beta}{\beta\gamma}} du\\ &=&[\frac{\beta\gamma}{\alpha\gamma+\alpha\beta+\beta\gamma}u^\frac{\alpha\gamma+\alpha\beta+\beta\gamma}{\beta\gamma}]_0^1\\ &=&\frac{\beta\gamma}{\alpha\gamma+\alpha\beta+\beta\gamma} \end{eqnarray*} }$
より、
${\displaystyle P(X=max(X,Y,Z))=P(X>max(Y,Z))=\frac{\beta\gamma}{\alpha\gamma+\alpha\beta+\beta\gamma}=\frac{1/\alpha}{1/\alpha+1/\beta+1/\gamma} }$
とわかります。

(3)

(2)が解ければボーナス問題です。
$U_1=u_1$ の下での条件付き確率 $P(X_1=max(X_1,...,X_n)|U_1=u_1)$ は、(2)での議論から、
「 $u_1^{\alpha_1}$ が、 $u_1^{\alpha_1}>u_2^{\alpha_2}$ ,..., $u_1^{\alpha_1}>u_n^{\alpha_n}$ を同時に満たす確率」と同じ意味なので、
(2)と同様に
$u_1^{\alpha_1/\alpha_2}$ ,..., $u_1^{\alpha_1/\alpha_n}$
を掛けあわせて
$P(X_1=max(X_1,...,X_n)|U_1=u_1)=u_1^{\alpha_1/\alpha_2+...+\alpha_1/\alpha_n}$
これを積分して
$P(X_1=max(X_1,...,X_n))=\int_0^1 P(X_1=max(X_1,...,X_n)|U_1=u_1)du_1=\frac{1/\alpha_1}{1/\alpha_1+...+1/\alpha_n}$
と求まります。

以下公式参考書等の紹介です(踏んで買って頂けたら嬉しいですが、そこそこお高いと思うので慎重にご判断下さい)

2012年度と2013年度については公式解説集があります。まだ施行歴が浅いため統計検定の過去問は2年分だけで、あとはRSS/JSSとの抱き合わせになっていますのでご注意を。

公式の参考書もあります。コンパクトに試験範囲が網羅されている反面、残念ながら定理の証明がかなり略されているのでもし使うなら試験範囲の確認用として、細かいところは他の本(統計学入門とか)で詰めた方が良いような気がします。僕は証明をしながら進めようとしたら全然読み終わらなくてヤバイヤバイ言ってます。

2015-09-15

PRML悪戦苦闘(演習問題 1.17)

PRML 問題演習計算問題ガンマ関数

いつまたやる気をなくすかわかりませんし、扱った問題で出てきた他章の問題については順番に依らず解こうと思います。

1.17

問題文

ガンマ関数は

$\Gamma(x)\equiv \int_0^\infty u^{x-1}e^{-u}du$

で定義される。

(1)部分積分を使って関係式 $\Gamma(x+1)=x\Gamma(x)$ を証明せよ。

(2) $\Gamma(1)=1$ を示し、xが整数なら $\Gamma(x+1)=x!$ となることを示せ。

解答

(1)
指示通りに部分積分します。
${\displaystyle \begin{eqnarray*} \Gamma(x+1) &=&\int_0^\infty u^{x}e^{-u}du\\ &=&\int_0^\infty u^{x}(-e^{-u})'du\\ &=&[ u^{x}(-e^{-u}) ]_0^\infty-\int_0^\infty (u^{x})'(-e^{-u})du\\ &=&\left[ -\frac{u^{x}}{e^u} \right]_0^\infty-\int_0^\infty (u^{x})'(-e^{-u})du\\ &=&(0-0)-\int_0^\infty x(u^{x-1})(-e^{-u})du\\ &=&x\int_0^\infty(u^{x-1})(-e^{-u})du=x\Gamma(x) \end{eqnarray*} }$
ということで示せました。

(2)
xが整数なら $\Gamma(x+1)=x!$ (*)となることを数学的帰納法で示します。

(i) $x=1$ のとき
${\displaystyle \begin{eqnarray*} \Gamma(1) &=&\int_0^\infty u^{0}e^{-u}du\\ &=&\int_0^\infty e^{-u}du\\ &=&\left[ -\frac{1}{e^u} \right]_0^\infty\\ &=&-(0-1)=1 \end{eqnarray*} }$

(ii)kを整数として、 $x=k+1$ のとき(*)が成立すると仮定し、 $x=k+2$ の場合に成り立つことを示します。

(1)を利用して、(*)の左辺を変形すると
${\displaystyle \begin{eqnarray*} \Gamma(k+2)=(k+2)\Gamma(k+1) \end{eqnarray*} }$

仮定より、
${\displaystyle \begin{eqnarray*} (k+2)\Gamma(k+1) &=&(k+2)(k+1)!=(k+2)! \end{eqnarray*} }$
これは(*)の右辺と一致します。
したがって(i),(ii)より、xが整数のとき、 $\Gamma(x+1)=x!$ が成り立ちます。

今回は易しめですね

2015-09-15

PRML悪戦苦闘(演習問題 2.10)

(追記)ガンマ関数の性質の証明を追加しましたblue0620.hatenablog.com

2.9が飛んでる気がするって？気のせいでしょう

2.10

問題文(適当に意訳)

ガンマ分布の性質 $\Gamma(x+1)=x \Gamma(x)$ を用いて、

$Dir({\bf \mu}|{\bf \alpha})=\frac{\Gamma(\alpha_0)}{\Gamma(\alpha_1)...\Gamma(\alpha_K)}\Pi_{k=1}^K \mu_k^{\alpha_k -1}$ (但し、 $\alpha_o=\sum_{k=1}^K \alpha_k$ )

のディリクレ分布の
(1)平均
${\bf E}[\mu_j]=\frac{\alpha_j}{\alpha_0}$
(2)分散
$var[\mu_j]=\frac{\alpha_j(\alpha_0-\alpha_j)}{\alpha_0 ^2(\alpha_0+1) }$
(3)共分散
$cov[\mu_j \mu_l]=-\frac{\alpha_j \alpha_l}{\alpha_0^2(\alpha_0+1)}$
の結果を導出せよ

方針

問題2.9で、K変数のディリクレ分布が正規化されていることが示されたので(示してないけど)、

$1\le j\le K$ なる整数 $j$ について
$\int_0 ^1 Dir({\bf \mu}|{\bf \alpha})d\mu_j=1$
が成り立ちます。
つまり
$\int_0 ^1 Dir({\bf \mu}|{\bf \alpha})d\mu_j=\frac{\Gamma(\alpha_0)}{\Gamma(\alpha_1)...\Gamma(\alpha_K)} \Pi_{k=1}^K \mu_k^{\alpha_k -1}$
より、
$\int_0 ^1 \Pi_{k=1}^K \mu_k^{\alpha_k -1}d\mu_j=\frac{\Gamma(\alpha_1)...\Gamma(\alpha_K)}{\Gamma(\alpha_0)}$

この正規化条件をうまく使うのがポイントみたいです。

解答

(1)
${\displaystyle \begin{eqnarray*} {\bf E}[\mu_j] &=&\int_0^1 \mu_j Dir(\bf{\mu}|\bf{\alpha})d\mu_j \\ &=&\int_0^1 \mu_j \frac{\Gamma(\alpha_0)}{\Gamma(\alpha_1)...\Gamma(\alpha_K)} \Pi_{k=1}^K \mu_k^{\alpha_k -1} d\mu_j\\ &=&\frac{\Gamma(\alpha_0)}{\Gamma(\alpha_1)...\Gamma(\alpha_K)} \color{red}{\int_0^1 \mu_j \Pi_{k=1}^K \mu_k^{\alpha_k -1} d\mu_j} \end{eqnarray*} }$

赤で表した部分について、
${\displaystyle \begin{eqnarray*} \int_0^1 \mu_j \Pi_{k=1}^K \mu_k^{\alpha_k -1} d\mu_j &=&\int_0^1 \mu_j(\mu_1^{\alpha_1 -1} \mu_2^{\alpha_2 -1}...\mu_K^{\alpha_K -1})d\mu_j\\ &=&\int_0^1 \mu_1^{\alpha_1 -1} \mu_2^{\alpha_2 -1}...\mu_j^{\alpha_j}...\mu_K^{\alpha_K -1} d\mu_j \end{eqnarray*} }$

したがって
$\int_0^1 \Pi_{k=1}^K \mu_k^{\alpha_k -1} d\mu_j$ において、 $\mu_j^{\alpha_j -1}$ を $\mu_j^{\alpha_j}$ に置き換えたものになります。

これを正規化項の関係式の右辺
${\displaystyle \begin{eqnarray*} \frac{\Gamma(\alpha_1)...\Gamma(\alpha_j)...\Gamma(\alpha_K)}{\Gamma(\alpha_0)}=\frac{\Gamma(\alpha_1)...\Gamma(\alpha_K)}{\Gamma(\alpha_1+\alpha_2+...+\alpha_K)} \end{eqnarray*} }$

についても同様に $\alpha_j$ を $\alpha_j +1$ に置き換えると、

${\displaystyle \begin{eqnarray*} \frac{\Gamma(\alpha_1)...\Gamma(\alpha_j+1)...\Gamma(\alpha_K)}{\Gamma(\alpha_1+\alpha_2+...+(\alpha_j+1)+...+\alpha_K)} =\frac{\Gamma(\alpha_1)...\Gamma(\alpha_j+1)...\Gamma(\alpha_K)}{\Gamma(\alpha_0+1)} \end{eqnarray*} }$

となります。したがってさっきの赤い部分に代入して

${\displaystyle \begin{eqnarray*} {\bf E}[\mu_j] &=&\frac{\Gamma(\alpha_0)}{\Gamma(\alpha_1)...\Gamma(\alpha_K)} \int_0^1 \mu_j \Pi_{k=1}^K \mu_k^{\alpha_k -1} d\mu_j\\ &=&\frac{\Gamma(\alpha_0)}{\Gamma(\alpha_1)...\Gamma(\alpha_K)} \frac{\Gamma(\alpha_1)...\Gamma(\alpha_j+1)...\Gamma(\alpha_K)}{\Gamma(\alpha_0+1)}\\ &=&\frac{\Gamma(\alpha_0)}{\Gamma(\alpha_1)...\Gamma(\alpha_K)} \frac{\Gamma(\alpha_1)...\color{red}{\alpha_j \Gamma(\alpha_j)}...\Gamma(\alpha_K)}{\color{red}{\alpha_0 \Gamma(\alpha_0)}}\\ &=&\frac{\alpha_j}{\alpha_0} \end{eqnarray*} }$
最後の赤字の変形が今回与えられたガンマ分布の性質 $\Gamma(x+1)=x \Gamma(x)$ を用いた箇所です。最後にきれいに求まると気持ち良いですね。

(2)
$var[\mu_j]={\bf E}[\mu_j^2]-{\bf E}[\mu_j]^2$ なので、いつもどおり ${\bf E}[\mu_j^2]$ を求めていきます。
式変形は(1)と同じ感じでやればよいですね。

${\displaystyle \begin{eqnarray*} {\bf E}[\mu_j^2] &=&\int_0^1 \mu_j^2 Dir(\bf{\mu}|\bf{\alpha})d\mu_j \\ &=&\int_0^1 \mu_j^2 \frac{\Gamma(\alpha_0)}{\Gamma(\alpha_1)...\Gamma(\alpha_K)} \Pi_{k=1}^K \mu_k^{\alpha_k -1} d\mu_j\\ &=&\frac{\Gamma(\alpha_0)}{\Gamma(\alpha_1)...\Gamma(\alpha_K)} \color{red}{\int_0^1 \mu_j^2 \Pi_{k=1}^K \mu_k^{\alpha_k -1} d\mu_j} \end{eqnarray*} }$

赤で表した部分について、

${\displaystyle \begin{eqnarray*} \int_0^1 \mu_j^2 \Pi_{k=1}^K \mu_k^{\alpha_k -1} d\mu_j &=&\int_0^1 \mu_j^2(\mu_1^{\alpha_1 -1} \mu_2^{\alpha_2 -1}...\mu_K^{\alpha_K -1})d\mu_j\\ &=&\int_0^1 \mu_1^{\alpha_1 -1} \mu_2^{\alpha_2 -1}...\mu_j^{\alpha_j+1}...\mu_K^{\alpha_K -1} d\mu_j \end{eqnarray*} }$

したがって
$\int_0^1 \Pi_{k=1}^K \mu_k^{\alpha_k -1} d\mu_j$ において、 $\mu_j^{\alpha_j -1}$ を $\mu_j^{\alpha_j+1}$ に置き換えたものになります。

これを正規化項の関係式の右辺についても同様に $\alpha_j$ を $\alpha_j +2$ に置き換えると、

${\displaystyle \begin{eqnarray*} \frac{\Gamma(\alpha_1)...\Gamma(\alpha_j+2)...\Gamma(\alpha_K)}{\Gamma(\alpha_1+\alpha_2+...+(\alpha_j+2)+...+\alpha_K)} =\frac{\Gamma(\alpha_1)...\Gamma(\alpha_j+2)...\Gamma(\alpha_K)}{\Gamma(\alpha_0+2)} \end{eqnarray*} }$

となります。したがって、さっきの赤い部分に代入して

${\displaystyle \begin{eqnarray*} {\bf E}[\mu_j^2] &=&\frac{\Gamma(\alpha_0)}{\Gamma(\alpha_1)...\Gamma(\alpha_K)} \int_0^1 \mu_j^2 \Pi_{k=1}^K \mu_k^{\alpha_k -1} d\mu_j\\ &=&\frac{\Gamma(\alpha_0)}{\Gamma(\alpha_1)...\Gamma(\alpha_K)} \frac{\Gamma(\alpha_1)...\Gamma(\alpha_j+2)...\Gamma(\alpha_K)}{\Gamma(\alpha_0+2)}\\ &=&\frac{\Gamma(\alpha_0)}{\Gamma(\alpha_1)...\Gamma(\alpha_K)} \frac{\Gamma(\alpha_1)...\color{red}{(\alpha_j+1) \alpha_j \Gamma(\alpha_j)}...\Gamma(\alpha_K)}{\color{red}{(\alpha_0+1)\alpha_0 \Gamma(\alpha_0)}}\\ &=&\frac{(\alpha_j+1)\alpha_j}{(\alpha_0+1)\alpha_0} \end{eqnarray*} }$
最後の式変形はガンマ関数の性質を2回繰り返して使っています。

あとは(1)で求めたものと併せて代入するだけです。
${\displaystyle \begin{eqnarray*} var[\mu_j] &=&{\bf E}[\mu_j^2]-{\bf E}[\mu_j]^2\\ &=&\frac{(\alpha_j+1)\alpha_j}{(\alpha_0+1)\alpha_0}-(\frac{\alpha_j}{\alpha_0})^2\\ &=&\frac{(\alpha_j+1)\alpha_j\alpha_0}{\alpha_0^2(\alpha_0+1)}-\left(\frac{\alpha_j^2(\alpha_0+1)}{\alpha_0^2(\alpha_0+1)}\right)^2 &=&\frac{\alpha_j(\alpha_0-\alpha_j)}{\alpha_0^2(\alpha_0+1)} \end{eqnarray*} }$

やったね！

(3)
$cov[\mu_j \mu_l]={\bf E}[\mu_j \mu_l]-{\bf E}[\mu_j]{\bf E}[\mu_l]$
を利用すると、 ${\bf E}[\mu_j \mu_l]$ を式変形すれば良いので、

${\displaystyle \begin{eqnarray*} {\bf E}[\mu_j \mu_l] &=&\int_0^1 \mu_j \mu_l Dir(\bf{\mu}|\bf{\alpha})d\mu_j \\ &=&\int_0^1 \mu_j \mu_l \frac{\Gamma(\alpha_0)}{\Gamma(\alpha_1)...\Gamma(\alpha_K)} \Pi_{k=1}^K \mu_k^{\alpha_k -1} d\mu_j\\ &=&\frac{\Gamma(\alpha_0)}{\Gamma(\alpha_1)...\Gamma(\alpha_K)} \color{red}{\int_0^1 \mu_j \mu_l \Pi_{k=1}^K \mu_k^{\alpha_k -1} d\mu_j} \end{eqnarray*} }$

赤で表した部分について、

${\displaystyle \begin{eqnarray*} \int_0^1 \mu_j \mu_l \Pi_{k=1}^K \mu_k^{\alpha_k -1} d\mu_j &=&\int_0^1 \mu_j \mu_l(\mu_1^{\alpha_1 -1} \mu_2^{\alpha_2 -1}...\mu_K^{\alpha_K -1})d\mu_j\\ &=&\int_0^1 \mu_1^{\alpha_1 -1} \mu_2^{\alpha_2 -1}...\mu_j^{\alpha_j+1}...\mu_l^{\alpha_l+1}...\mu_K^{\alpha_K -1} d\mu_j \end{eqnarray*} }$

これを正規化項の関係式の右辺についても同様に $\alpha_j$ を $\alpha_j +1$ に、 $\alpha_l$ を $\alpha_l +1$ にそれぞれ置き換えると、

${\displaystyle \begin{eqnarray*} \frac{\Gamma(\alpha_1)...\Gamma(\alpha_j+1)...\Gamma(\alpha_l+1)...\Gamma(\alpha_K)}{\Gamma(\alpha_1+\alpha_2+...+(\alpha_j+1)+...+(\alpha_l+1)+...+\alpha_K)} =\frac{\Gamma(\alpha_1)...\Gamma(\alpha_j+1)...\Gamma(\alpha_l+1)...\Gamma(\alpha_K)}{\Gamma(\alpha_0+2)} \end{eqnarray*} }$

となります。したがって、さっきの赤い部分に代入して

${\displaystyle \begin{eqnarray*} {\bf E}[\mu_j \mu_l] &=&\frac{\Gamma(\alpha_0)}{\Gamma(\alpha_1)...\Gamma(\alpha_K)} \int_0^1 \mu_j \mu_l \Pi_{k=1}^K \mu_k^{\alpha_k -1} d\mu_j\\ &=&\frac{\Gamma(\alpha_0)}{\Gamma(\alpha_1)...\Gamma(\alpha_K)} \frac{\Gamma(\alpha_1)...\Gamma(\alpha_j+1)...\Gamma(\alpha_l+1)...\Gamma(\alpha_K)}{\Gamma(\alpha_0+2)}\\ &=&\frac{\Gamma(\alpha_0)}{\Gamma(\alpha_1)...\Gamma(\alpha_K)} \frac{\Gamma(\alpha_1)...\color{red}{\alpha_j\Gamma(\alpha_j)\alpha_l\Gamma(\alpha_l)}...\Gamma(\alpha_K)}{\color{red}{(\alpha_0+1)\alpha_0 \Gamma(\alpha_0)}}\\ &=&\frac{\alpha_j\alpha_l}{(\alpha_0+1)\alpha_0} \end{eqnarray*} }$
やったね！

完全に式変形がワンパターン化してますね。次はせっかくなので少し戻って今回出てきたガンマ関数の性質の証明(演習1.17)を見ていこうと思います。

2015-09-15

PRML悪戦苦闘(演習問題 2.8)

PRML 計算問題加法定理

僕の怠惰(と学会準備)のために随分間が空いてしまいました。

PRML勉強会で、文中に出てくる演習問題はまじめに取り組んでいたので懐かしかった。このブログは頑張って続けてほしい。｜PRML悪戦苦闘(演習問題 2.7) - blueblog URL

2015-09-08 20:55:10 via YoruFukurou

大変ありがたいことに励まして下さった方がいらしたので、また少しずつ再開したいと思います。

2.8

問題文

同時確率が $p(x,y)$ であるような2つの変数 $x,y$ を考える。これについて、次の2つの結果を証明せよ。

(1)
${ \displaystyle \begin{eqnarray*} {\bf E}[x]={\bf E}_y [{\bf E}_x [x|y]] \end{eqnarray*} }$

(2)
${ \displaystyle \begin{eqnarray*} var[x]={\bf E}_y [var_x[x|y]]+var_y[{\bf E}_x [x|y]] \end{eqnarray*} }$

解答

(1)
${ \displaystyle \begin{eqnarray*} {\bf E}_x [x|y]&=&\int p(x|y)x dx\\ &=&\int \frac{p(x,y)}{p(y)}x dx \end{eqnarray*} }$

これは単に「yが起こる事象下で、xが起こる確率」すなわちベイズの条件付き確率の定義を2通りの形で表しただけですね。
最後の式の分母がxを含まないことを下の変形で利用します。

${ \displaystyle \begin{eqnarray*} {\bf E}_y [{\bf E}_x [x|y]] &=&\int p(y) \int p(x|y)x dx dy\\ &=&\int p(y) \int \frac{p(x,y)}{p(y)}x dx dy\\ &=&\int \int p(x,y)x dx dy\\ &=&\color{red}{\int x\int p(x,y)dy dx}\\ &=&\color{red}{\int x p(x)dx}={\bf E}[x] \end{eqnarray*} }$

ここで赤色の式変形に注意しましょう(基本なのですが一応)。これは確率の加法定理を表します。
離散値の場合これは $p(x)=\sum_y p(x,y)$ で表わせ、xとyの同時確率を全てのyについて足し合わせれば、xが起こる確率になることを意味します。
(大きなサイコロxと小さなサイコロyを同時に振るとき、xが1となる確率は
「x=1かつy=1」,「x=1かつy=2」,「x=1かつy=3」…,「x=1かつy=6」となる確率の和ですよね？)

(2)
${ \displaystyle \begin{eqnarray*} {\bf E}_y[var_x[x|y]]+var_y[{\bf E}_x[x|y]]&=&{\bf E}_y[ {\bf E}_x[(x|y)^2]-{\bf E}_x [x|y]^2]+{\bf E}_y[{\bf E}_x[x|y]^2]-{\bf E}_y [{\bf E}_x[x|y]]^2\\ &=&{\bf E}_y[{\bf E}_x[(x|y)^2]]-{\bf E}_y[{\bf E}_x [x|y]]^2\\ &=&\int p(y)\int x^2 \frac{p(x,y)}{p(y)}dx dy-{\bf E}_y[{\bf E}_x [x|y]]^2\\ &=&\color{red}{\int \int x^2 p(x,y)dx dy}-{\bf E}[x]^2 \\ &=&\color{red}{\int x^2 p(x)dx} - {\bf E}[x]^2\\ &=&{\bf E}[x^2]-{\bf E}[x]^2=var[x] \end{eqnarray*} }$

ここでの赤い式変形も加法定理です。

基本的なベイズの定理と加法定理を式変形でいかにさらっと使えるか、という問題だと思います。僕はそんなにさらっと書けませんでしたけど（泣）

2015-07-31

PRML悪戦苦闘(演習問題 2.7)

2.7

問題文(適当に要約)

$\mu$ の事前分布がベータ分布
$Beta(\mu|a,b)=\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \mu^{a-1} (1-\mu)^{b-1}$
である二項分布
$\sum_{m=0}^N \begin{pmatrix}N\\m\end{pmatrix} \mu^m(1-\mu)^{N-m}$
にしたがう確率変数 $x$ を考える。

ここで $x=1$ の事象が $m$ 回、 $x=0$ の事象が $l$ 回生じたとする。このとき $\mu$ の事後平均が事前平均と $\mu$ の最尤推定量の間の値になることを示せ。

方針

事前分布の平均は、ベータ分布の超パラメータ $a,b$ を用いて
${ \displaystyle \frac{a}{a+b} }$
と表せます。
二項分布を尤度関数として、事前分布にベータ分布を選択すると、事後分布もまたベータ分布となることをPRMLの本編2.1.1で学びました。
これは共役性と呼ばれる重要な性質で、事後分布の超パラメータは $a$ を $a+m$ に、 $b$ を $b+l$ に置き換えたものになります。

従って事後平均は、
${ \displaystyle p(x=1|D)=\frac{a+m}{a+m+b+l} }$
と表せます。

また、データからの最尤推定量は単に観察値に基づいた確率を求めればよく、 $m+l$ 回の試行の中で $x=1$ の事象が起こる確率
${ \displaystyle \frac{m}{m+l} }$
となります。

事後平均が事前平均と最尤推定量の間の値になることは、 $0\leq \lambda\leq 1$ として
${ \displaystyle \frac{a+m}{a+m+b+l}=\lambda \frac{a}{a+b}+(1-\lambda) \frac{m}{m+l} }$
を示せば良いですね(ここまではもとの問題文に書いてあります)。

解答

${ \displaystyle \begin{eqnarray*} \lambda \frac{a}{a+b}+(1-\lambda) \frac{m}{m+l}&=& \frac{a\lambda}{a+b}+ \frac{m(1-\lambda)}{m+l}\\ &=&\frac{a\lambda(m+l)+m(1-\lambda)(a+b)}{(a+b)(m+l)}\\ &=&\frac{al\lambda + m(a+b-b\lambda)}{(a+b)(m+l)}\\ &=&\frac{(al-bm)\lambda + m(a+b)}{(a+b)(m+l)} \end{eqnarray*} }$

したがって、
${ \displaystyle \begin{eqnarray*} \frac{m+a}{a+b+m+l}&=&\frac{(al-bm)\lambda + m(a+b)}{(a+b)(m+l)}\\ \end{eqnarray*} }$

これを $\lambda$ に関して式変形して
${ \displaystyle \begin{eqnarray*} (al-bm)\lambda&=&\{\frac{m+a}{a+b+m+l}-\frac{m}{m+l}\}*(m+l)(a+b)\\ &=&\frac{al-bm}{(a+b+m+l)(m+l)}(m+l)(a+b)\\ \end{eqnarray*} }$

したがって、
${ \displaystyle \begin{eqnarray*} \lambda=\frac{m+a}{a+b+m+l}\\ \end{eqnarray*} }$

a,b,m,lはいずれも非負なので $0\leq \frac{m+a}{a+b+m+l}\leq 1$ が示せます。

また、観察値m,lが大きくなると事後平均が最尤推定値に近づいていくことも理解できます。

2015-07-29

PRML悪戦苦闘(演習問題 2.6)

2.6

問題文

ベータ分布の平均・分散・およびモードがそれぞれ

(1) ${\bf E}[\mu]=\frac{a}{a+b}$

(2) $var[\mu]=\frac{ab}{(a+b)^2 (a+b+1)}$

(3) $mode[\mu]=\frac{a-1}{a+b-2}$

になることを示せ

方針

(1),(2)
代入してひたすら変形。xが自然数のときのガンマ関数の性質 $\Gamma(x+1)=x!$ を使える形にしていきます。

(3)
最頻値なので、ベータ分布の式を $\mu$ について微分して、イコール0になる条件のうち、極大値を取る時が答えです。

解答

(1)
${ \displaystyle \begin{eqnarray*} {\bf E}[\mu]&=&\int_0^1 \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}\mu^{a-1}(1-\mu)^{b-1}\mu d\mu\\ &=&\int_0^1 \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}\mu^{a}(1-\mu)^{b-1} d\mu\\ &=&\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}*\frac{\Gamma(a+1)\Gamma(b)}{\Gamma(a+b+1)}\\ &=&\frac{\Gamma(a+b)\Gamma(a+1)}{\Gamma(a)\Gamma(a+b+1)}\\ &=&\frac{(a+b+1)!a!}{(a-1)!(a+b)!}\\ &=&\frac{a!}{(a-1)!(a+b)}=\frac{a}{(a+b)} \end{eqnarray*} }$

(2)
${ \displaystyle \begin{eqnarray*} {\bf E}[\mu^2]&=&\int_0^1 \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}\mu^{a-1}(1-\mu)^{b-1}\mu^2 d\mu\\ &=&\int_0^1 \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}\mu^{a+1}(1-\mu)^{b-1} d\mu\\ &=&\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}*\frac{\Gamma(a+2)\Gamma(b)}{\Gamma(a+b+2)}\\ &=&\frac{(a+b+1)!(a+1)!}{(a-1)!(a+b+1)!}=\frac{(a+1)a}{(a+b)(a+b+1)} \end{eqnarray*} }$

${ \displaystyle \begin{eqnarray*} var[\mu]&=&{\bf E}[\mu^2]-{\bf E}[\mu]^2\\ &=&\frac{(a+1)a}{(a+b)(a+b+1)}-\frac{a^2}{(a+b)^2}\\ &=&\frac{a(a+1)(a+b)-a^2(a+b+1)}{(a+b)^2 (a+b+1)}\\ &=&\frac{ab}{(a+b)^2 (a+b+1)} \end{eqnarray*} }$

(3)

ベータ分布の式は $Beta(\mu|a,b)=\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \mu^{a-1} (1-\mu)^{b-1}$ ですが、 $\mu$ の変化による増減を考えるので $\mu^{a-1} (1-\mu)^{b-1}$ の部分のみに注目します。

${ \displaystyle \begin{eqnarray*} (\mu^{a-1} (1-\mu)^{b-1})'&=&(a-1)\mu^{a-2}(1-\mu)^{b-1}-\mu^{a-1}(b-1)(1-\mu)^{b-2}\\ &=&(1-\mu)^{b-2}\mu^{a-2}\{(a-1)(1-\mu)-\mu(b-1) \}\\ &=&(1-\mu)^{b-2}\mu^{a-2}\{a-1-\mu a-\mu b+2\mu \}\\ &=&(1-\mu)^{b-2}\mu^{a-2}\{a-1-\mu(2+a+b)\}=0 \end{eqnarray*} }$

求めたいのは極大値をとるとき(最頻値をとるとき)なので
${ \displaystyle mode[\mu]=\frac{a-1}{2+a+b} }$
と求まります。