2015年10月19日

魔方陣トーク

自宅での話。

次男坊が3×3の魔方陣を考えている様子。

と、そこへ入ってきた兄。6年生なので3x3の魔方陣くらいは簡単に答えてしまいます。やりたがりなので、それを見つけるやいなや弟から紙片を取り上げたところでストップをかけました。6年生にはその問題は簡単すぎるから、弟君にやらせてあげなさい。そのかわり次の問題を考えてみよう、と言って出したのが:

「3x3の魔方陣の中心の数字は5であり、それ以外にないことを証明せよ」(但し、3x3の魔方陣は通常のもの、すなわち各マスに1〜9の数字が一度ずつ入り、縦横斜め全ての三つ組みの和が同じであるようなものとした。)

5が中心に来ることは経験的にわかっているようですが、それ以外の数字は絶対にはいらないことを証明するという問題です。彼が出した答えは次。

魔方陣を次のように表すと、
\begin{array}{ccc}
i & b & c \\
h & a & d \\
g & f & e
\end{array}

1) まず、魔方陣は全て足して45(∵1+2+...+9=45)で、縦横斜めどれを足しても同じ数だから、三つ組みの和は15でなければならない。

2) 次に、aは偶数ではない。なぜなら、三つ組みの和が15で奇数なので、aが偶数だとaを通る縦横斜めの三つ組みの、a以外の数字は奇数と偶数が一つずつ入ることになる。すると、例えば(i,a,e), (h,a,d), (g,a,c)の三つに注目すると、(i,h,g)を足して15であるからこのうち一つまたは三つが奇数。ということは、(c,d,e)は一つまたは三つが偶数となるが、これは(c,d,e)を足して15即ち奇数となることと矛盾する。よって、aは奇数である。

3) aは奇数であるが、まず9ではない。なぜなら、もしa=9だとすると、b〜iのうちどれかに8が入るが、8と9の入る三つ組み(8,9,*)について8+9は既に15を超えているので、魔方陣の定義と矛盾する。次に、7ではない。これも同様に、b〜iのどこかに9が入るため。また、a=1だとすると、b〜iのうちどこかに2が入るが、その三つ組み(2,1,*)は残りのどの数字を入れても15には届かない。a=3も(1,3,*)について同様。

4) 以上より、aに入るのは奇数であり、また(9,7,3,1)のどれでもない。一方で、a=5である魔方陣は存在するので、この魔方陣の中心は5であり、それ以外ではない。

なるほど、背理法をうまく活用しました。実は2)無しで、3)の方法だけでも十分なのですが(a=6以上だと(9,a,*)の三つ組みが15を超え、a=4以下だと(1,a,*)の三つ組みが15に届かない)、試行錯誤の順番としては悪くないです。

そのあと、もうちょっと面白い方法はない?と言ったら出てきたのが次の方法。

\begin{array}{ccc}
G \searrow & D \downarrow & E \downarrow & F \downarrow & \swarrow H \\
A \rightarrow & i & b & c \\
B \rightarrow & h & a & d \\
C \rightarrow & g & f & e
\end{array}

魔方陣の各三つ組みにA〜Hの名前をつけました。A,B,Cが行、D,E,Fが列、G,Hが斜めです。これらはすべて、三つ数字を和して15です。すると、これらをうまく足し引きして、次の式が成立します。

\begin{align}
A+C+D+F-G-H+3 \times a=a+b+c+d+e+f+g+h+i
\end{align}

これは、A+Cで横、D+Fで縦の、真ん中を通らない行と列をたすと、真ん中以外のb〜iが一つずつと、角がもう一つずつたされて、角をたしすぎなので斜めを引いたら、今度は真ん中が2回引かれ過ぎなので3つたして、そうするとこれでa〜iの全部が1回ずつたされたという式です。ところで、a〜iまでたすと45、また、A〜Hはすべて15なので、A+C+D+F-G-Hは30。よって、

\begin{align}
30+3 \times a=45
\end{align}

となり、めでたくa=5が出てきました。


posted by jinya at 13:33| Comment(0) | 日記 | このブログの読者になる | 更新情報をチェックする

2015年10月10日

回帰モデルのギブスサンプラーを導出する(その2)

回帰モデルのギブスサンプラーを導出します。教科書は「マルコフ連鎖モンテカルロ法」(豊田秀樹著)。こちら http://jinya.seesaa.net/article/427453869.html の続きです。

回帰モデル
\begin{align}
y=X \beta + \varepsilon, \: \varepsilon \sim N(0,\tau^{-1}I_n)
\end{align}
但し、サンプル数は\(n\)、説明変数の数は\(K-1\)個(定数項を入れて\(K\)次元)とする。\(I_n\)は\(n\times n\)の単位行列。なお、教科書では\(n\)の所を\(I\)としていますが、単位行列の\(I\)と間違いやすいのでここでは\(n\)を使います。

このモデルに対して、事前分布を次のように採ります。
\begin{align}
[\beta | \tau,X,y] \sim N(b_0, \tau^{-1}B_0), \: [\tau,\beta,X,y] \sim Ga(n_0 / 2, n_0 S_0 /2)
\end{align}
この事前分布から事後分布を求めます。

【事前分布(同時)\(\pi(\beta,\tau)\)】
\begin{align}
&\pi(\beta,\tau ) \\
&\propto \tau^{K/2} \exp \left( -\frac{1}{2}(\beta-b_0)'(\tau^{-1}B_0)^{-1}(\beta-b_0) \right) \tau^{n_0/2-1}\exp \left( -\frac{n_0S_0}{2} \tau \right) \\
&=\tau^{(n_0+K)/2-1} \exp \left( -\frac{\tau}{2} \left( n_0S_0 + (\beta-b_0)'(\tau^{-1}B_0)^{-1}(\beta-b_0) \right) \right)
\end{align}

【尤度関数】
\begin{align}
\pi(y|\beta,\tau,X) \propto \tau^{n/2}\exp \left( -\frac{\tau}{2}(y-X\beta)'(y-X\beta) \right)
\end{align}

【事後分布】
事後分布は定数の差を除いて事前分布と尤度関数とのかけ算なので、
\begin{align}
\pi(&\beta,\tau | X,y) \propto \tau^{(n_0+n+K)/2-1} \\
& \times \exp \left( -\frac{\tau}{2} \left( \underline{ n_0S_0 + (\beta-b_0)'B_0^{-1}(\beta-b_0) + (y-X\beta)'(y-X\beta)}_{\fbox{1}} \right) \right)
\end{align}

\(\fbox{1}\)について、
\begin{align}
\fbox{1} = n_0S_0 + (\beta-b_0)'B_0^{-1}(\beta-b_0) + \underline{(y-X\beta)'(y-X\beta)}_{\fbox{2}}
\end{align}

\(\fbox{2}\)について、最尤推定された\(\beta\)を\(\hat{\beta}\)とすると、\( y = \hat{y} + \varepsilon, \: \hat{y} = X\hat{\beta} \) を用いて、
\begin{align}
& (y-X\beta)'(y-X\beta) \\
&= ( \varepsilon + \hat{y} - X\beta)'( \varepsilon + \hat{y} - X\beta) \\
&= \varepsilon' \varepsilon + (\hat{y} - X\beta)'(\hat{y} - X\beta) + \varepsilon'(\hat{y} - X\beta) + (\hat{y} - X\beta)' \varepsilon \\
&= S_e + (\hat{\beta}-\beta)'X'X(\hat{\beta}-\beta) + \underline{\varepsilon'X}(\hat{\beta}-\beta)+ (\hat{\beta}-\beta)'\underline{X'\varepsilon} \\
&= S_e + Q(\beta)
\end{align}
但し、\(S_e=\varepsilon'\varepsilon\), \(Q(\beta)=(\hat{\beta}-\beta)'X'X(\hat{\beta}-\beta)\)とした。下線部はゼロ。

戻って、
\begin{align}
\fbox{1}=n_0S_0 + S_e + \underline{(\beta-b_0)'B_0^{-1}(\beta-b_0) + Q(\beta)}_\fbox{3}
\end{align}

\(B_1^{-1}=B_0^{-1}+X'X\)及び\( b_1=B_1(B_0^{-1}b_0 + X'X\hat{\beta}) \)とすると、

\begin{align}
\fbox{3}
&=(\beta-b_0)'B_0^{-1}(\beta-b_0)+(\hat{\beta}-\beta)'X'X(\hat{\beta}-\beta) \\
&= \beta'B_0^{-1}\beta - \beta'B_0^{-1}b_0 - b_0'B_0^{-1}\beta + b_0'B_0^{-1}b_0 \\
&+ \hat{\beta}'X'X\hat{\beta} - \beta'X'X\hat{\beta} - \hat{\beta}'X'X\beta + \beta'X'X\beta \\
&= \beta'(B_0^{-1}+X'X)\beta - \beta'(B_0^{-1}b_0+X'X\hat{\beta}) \\
& - (b_0'B_0^{-1}+\hat{\beta}'X'X)\beta + \underline{b_0'B_0^{-1}b_0 + \hat{\beta}'X'X\hat{\beta}}_{\fbox{4}} \\
&= \beta'B_1^{-1}\beta - \beta'B_1^{-1}b_1 - b_1'B_1^{-1}\beta + \fbox{4}
\end{align}

ここで、\(B_0, X'X\)は対称行列なので\(B_1\)も対称行列であり、さらにそれらの逆行列も対称行列である性質を使っています。

\begin{align}
\fbox{4}&=b_0'B_0^{-1}b_0 + \hat{\beta}'X'X\hat{\beta} \\
&= b_0' B_1^{-1} B_1 B_0^{-1} b_0 + \hat{\beta}' B_1^{-1} B_1 X'X \hat{\beta} \\
&= b_0' (B_0^{-1}+X'X) B_1 B_0^{-1} b_0 + \hat{\beta}' (B_0^{-1}+X'X) B_1 X'X \hat{\beta} \\
&= b_0' B_0^{-1} B_1 B_0^{-1} b_0 + b_0' X'X B_1 B_0^{-1} b_0
+ \hat{\beta}' B_0^{-1} B_1 X'X \hat{\beta} + \hat{\beta}' X'X B_1 X'X \hat{\beta} \\
&= (b_0'B_0^{-1}+\hat{\beta}'X'X)B_1(B_0^{-1}b_0+X'X\hat{\beta}) \\
& -b_0'B_0^{-1}B_1X'X\hat{\beta} -\hat{\beta}' X'X B_1 B_0^{-1} b_0
+ b_0' X'X B_1 B_0^{-1} b_0 + \hat{\beta}' B_0^{-1} B_1 X'X \hat{\beta} \\
&= b_1'B_1^{-1}b_1 + (b_0'-\hat{\beta}')X'XB_1B_0^{-1}b_0 - (b_0'-\hat{\beta})B_0^{-1}B_1X'X\hat{\beta} \\
&= b_1'B_1^{-1}b_1 + (b_0-\hat{\beta})'B_0^{-1}B_1X'X(b_0-\hat{\beta})
\end{align}

よって、
\begin{align}
\fbox{3}&=\beta'B_1^{-1}\beta - \beta'B_1^{-1}b_1 - b_1'B_1^{-1}\beta + b_1'B_1^{-1}b_1 + (b_0-\hat{\beta})'B_0^{-1}B_1X'X(b_0-\hat{\beta}) \\
&= (\beta-b_1)'B_1^{-1}(\beta-b_1) + (b_0-\hat{\beta})'B_0^{-1}B_1X'X(b_0-\hat{\beta})
\end{align}

\begin{align}
\fbox{1}&= n_0S_0 + S_e + \fbox{3} \\
&= n_0S_0 + (\beta-b_1)'B_1^{-1}(\beta-b_1)
+ \underline{S_e + (b_0-\hat{\beta})'B_0^{-1}B_1X'X(b_0-\hat{\beta})}_\fbox{5}
\end{align}

\begin{align}
\fbox{5}&=S_e + (b_0-\hat{\beta})'B_0^{-1}B_1X'X(b_0-\hat{\beta}) \\
&= \underline{S_e + (\hat{\beta}-b_0)'B_0^{-1}B_1X'X \hat{\beta}}_\fbox{6}
+ \underline{(b_0-\hat{\beta})'B_0^{-1}B_1X'X b_0}_\fbox{7}
\end{align}

ここで、\(B_1^{-1}=B_0^{-1}+X'X \)より、\(B_0^{-1}B_1=(B_1^{-1}-X'X)B_1 = I-X'XB_1\)を用いて

\begin{align}
\fbox{6}&=S_e + (\hat{\beta}-b_0)'B_0^{-1}B_1X'X \hat{\beta} \\
&= S_e + \hat{\beta}'B_0^{-1}B_1X'X \hat{\beta} - b_0'B_0^{-1}B_1X'X \hat{\beta} \\
&= S_e + \hat{\beta}' (I-X'XB_1) X'X \hat{\beta} - b_0'B_0^{-1}B_1X'X \hat{\beta} \\
&= S_e + \hat{\beta}'X'X\hat{\beta} - \hat{\beta}' X'XB_1 X'X \hat{\beta} - b_0'B_0^{-1}B_1X'X \hat{\beta} \\
&= S_e + \hat{\beta}'X'X\hat{\beta} - (\hat{\beta}' X'X + b_0'B_0^{-1})B_1X'X \hat{\beta} \\
&= S_e + \hat{\beta}'X'X\hat{\beta} - (B_1(B_0^{-1}b_0 + X'X\hat{\beta}))' X'X \hat{\beta} \\
&= S_e + \hat{\beta}'X'X\hat{\beta} - b_1' X'X \hat{\beta} \\
&= S_e + (\hat{\beta}-b_1)'X'X\hat{\beta} \\
&= \varepsilon'\varepsilon + (\hat{\beta}-b_1)'X'X\hat{\beta} \\
&= (\varepsilon' + (\hat{\beta}-b_1)'X')(\varepsilon + X\hat{\beta}) \\
&= (\varepsilon + \hat{y} -X b_1)'(\varepsilon + \hat{y}) \\
&= (y -X b_1)'y
\end{align}

\begin{align}
\fbox{7}&=(b_0-\hat{\beta})'B_0^{-1}B_1X'X b_0 \\
&= b_0'B_0^{-1}B_1X'X b_0 -\hat{\beta}'B_0^{-1}B_1X'X b_0 \\
&= b_0'B_0^{-1}B_1(B_1^{-1}-B_0^{-1}) b_0 -\hat{\beta}'X'XB_1B_0^{-1} b_0 \\
&= b_0'B_0^{-1}b_0 - b_0'B_0^{-1}B_1 B_0^{-1} b_0 -\hat{\beta}'X'XB_1B_0^{-1} b_0 \\
&= b_0'B_0^{-1}b_0 - (b_0'B_0^{-1} + \hat{\beta}'X'X) B_1 B_0^{-1} b_0 \\
&= b_0'B_0^{-1} b_0 - b_1' B_0^{-1} b_0 \\
&= (b_0-b_1)'B_0^{-1} b_0 \\
\end{align}

よって、\(\fbox{1}\)まで戻って、また、\(n_1=n_0+n\)、\(n_1S_1=n_0S_0 + (y -X b_1)'y + (b_0-b_1)'B_0^{-1} b_0 \)とすると、

\begin{align}
\fbox{1} &= n_0S_0 + (\beta-b_1)'B_1^{-1}(\beta-b_1) + \fbox{6} + \fbox{7} \\
&= \underline{n_0S_0 + (y -X b_1)'y + (b_0-b_1)'B_0^{-1} b_0} + (\beta-b_1)'B_1^{-1}(\beta-b_1) \\
&= n_1S_1 + (\beta-b_1)'B_1^{-1}(\beta-b_1)
\end{align}

さて、事後分布はどうだったかというと、
\begin{align}
&\pi(\beta,\tau | X,y)
\propto \tau^{(n_0+n+K)/2-1} \exp \left( -\frac{\tau}{2} \left( \fbox{1} \right) \right) \\
&= \underline{\tau^{K/2} \exp \left( -\frac{\tau}{2}(\beta-b_1)'B_1^{-1}(\beta-b_1) \right)}_\fbox{8}
\times \underline{\tau^{n_1/2-1} \exp \left( -\frac{\tau}{2} n_1S_1 \right)}_\fbox{9}
\end{align}

すると、\(\fbox{8}\)は正規分布\([\beta | \tau,X,y] \sim N(b_1,\tau B_1)\)の密度関数、\(\fbox{9}\)はガンマ分布\( [\tau | \beta,X,y] \sim Ga(n_1/2,n_1S_1/2)\)の密度関数であり、事後分布を解析的に書き下すことができた。さらに、\(\tau\)の従うガンマ分布の母数には\(\beta\)が入っていないので、\(\tau\)のサンプリングには\(\beta\)を必要としないこともわかる。
posted by jinya at 20:35| Comment(0) | 日記 | このブログの読者になる | 更新情報をチェックする

2015年10月08日

回帰モデルのギブスサンプラーを導出する(その1)

MCMCの教科書「マルコフ連鎖モンテカルロ法」(豊田秀樹著)はMCMCの中で何が起っているのかを腑に落とすのにわかりやすい教科書ですが、細かいデリベーションについては端折ってある場所もあって、そのなかで一カ所、回帰モデルのギブスサンプラー導出の記述で気になったところがあったので、詳細な式の導出をしてみます。場所は7.2節、「回帰モデルにおけるベイズ分析」です。

ところで、共役事前分布としてよくある形式は、BUGSっぽく書くと(但し、簡単のためモデル例は一次元)

model
{
for (i in 1:N){
Y[i] ~ dnorm(beta * X[i], tau);
}
tau ~ dgamma(0.001,0.001);
beta ~ dnorm(0,0.001);
}

これは回帰モデル
\begin{align}
Y=X \beta + \varepsilon, \: \varepsilon \sim N(0,1/\tau)
\end{align}
に対して、事前分布を
\begin{align}
&\beta \sim N(\mu_0, 1/\tau_0), \: \tau \sim Ga(s_0, r_0) \\
&(\mu_0=0, \tau_0=0.001, s_0=0.001, r_0=0.001)
\end{align}
で定義したモデルになります。

このとき、事後分布は解析的に求まって(共役事前分布を使うので)、
\begin{align}
\beta | \tau,X,Y \sim N(\mu_1, 1/\tau_1), \: \tau | \beta,X,Y \sim Ga(s_1,r_1)
\end{align}
但し(サンプル数を\(n\)とし、また簡単のため一次元)
\begin{align}
& \tau_1 = \tau_0 + \tau X'X, \: \mu_1 = (\tau_0 \mu_0 + \tau X'Y)/\tau_1 \\
& s_1 = s_0 + n/2, r_1 = r_0 + (Y-X\beta)'(Y-X\beta)/2
\end{align}

ギブスサンプラーですから、当然\(\beta\)の事後分布の中には\(\tau\)があり、また\(\tau\)の事後分布の中には\(\beta\)があります。

ところで、上述の教科書のモデルとはこれとはちょっと違っていて、事前分布が(一次元では)
\begin{align}
\beta | \tau \sim N(\mu_0, 1/(\tau_0 \tau)), \: \tau \sim Ga(s_0, r_0)
\end{align}
という形になっています。(変数は適宜読み替えてください。また、あとで教科書通りの変数名で再度、多次元版のデリベーションをします。)つまり、\(\beta\)の事前分布が既に\(\tau\)に依存しています。BUGSっぽいモデルで書くと、

model
{
for (i in 1:N){
Y[i] ~ dnorm(beta * X[i], tau);
}
tau ~ dgamma(0.001,0.001);
beta ~ dnorm(0,tau*0.001);
}

となって、事前分布のdnormの中にtauが入っています。これら二つのモデルは、事前分布が異なるので推定結果も異なります。しかし、そもそも事前分布に共役事前分布を使っている時点で事後分布は歪んでいるわけですから、そこはあまり気にしないことにします。

この修正された事前分布を用いると何が起るかというと、\(\tau\)の事後分布が\(\beta\)に依存しなくなります。もちろん観測値\((X,Y)\)には依存するので心配はいらないのですが、\(\tau\)が\(\beta\)に依存しないと言うことは、単純にギブスサンプリングでの計算量が半分ほどに減るわけで、そもそも計算量を減らすためにメトロポリス-ヘイスティングスではなくギブスサンプラーを使っていることを考えると、効率がより向上することになります。

長くなってきたので、続きは次回。

posted by jinya at 13:08| Comment(0) | 日記 | このブログの読者になる | 更新情報をチェックする
×

この広告は180日以上新しい記事の投稿がないブログに表示されております。