数式の夢とコンピュータの現実

bashにてteeとプロセス置換を使ってパイプを分岐する方法の備忘録です。

bashのパイプラインはテキスト処理に便利ですね。ビッグデータが普通に取り扱われるようになった昨今ですが、システムに組み込み済みの処理はともかく、アドホックなデータ分析業務の場合はまだまだテキストファイルでデータを受け渡しすることが多いです。そうすると、データが巨大な場合はそれをちょっと確認するだけでも一苦労もふた苦労もあったりして、そんなときにbashのテキスト処理系コマンドが威力を発揮します。

ここで、途中まで同じで、後ろが異なる処理を複数やりたい場合がありました。例えば、

zcat file0.gz | ./script0 | ./script1 > file1
zcat file0.gz | ./script0 | ./script2 > file2
zcat file0.gz | ./script0 | ./script3 > file3

のようなケースです。

ここで、file0.gzが小さくて、script0が一瞬で終わるようなものだったらどうってことないのですが、数ギガとか数十ギガという巨大なファイルを取り扱う場合には、このscript0を三回繰り返すのが勿体ない。かといって、script0を適用した結果を保存しておくにはディスクスペースが勿体ない。

そこで、teeを使ってパイプラインを分岐し、それぞれ結果を出力します。

zcat file0.gz | ./script0 | tee >(./script1 > file1) \
    | tee >(./script2 | file2) | ./script3 > file3

こうすると、file0の展開ファイルを作る必要が無く、また、script0を適用する回数も一回で済み、やりたい三つの作業を実行することができます。繰り返しになりますが、file0.gzが小さかったり、script0の処理が軽かったりすればなにもこんな面倒なことをしなくて済むのですが、巨大で重い場合はこれだけでかなりの時間節約になります。

なお、各scriptがメモリを大きく食う場合はメモリ不足に要注意です。そもそもパイプラインは一行単位のテキスト処理であるケースが多いですから、その場合はあまり気にしなくても大丈夫だろうと思います。


ちなみにこれを何に使ったかというと、あるPOSデータを記した巨大なcsvファイルの整形です。このファイルには、商品情報と店舗情報がマスターになっておらず、一つのテーブルにまとめて書いてありましたので、

zcat pos.csv.gz | nkf -Sw | sed -e 's/,/\t/g' -e 's/\"//g' \
    | tee >(cut -f 3,4,5,6 | sort -u > item.txt) \
    | tee >(cut -f 7,8,9,10 | sort -u > shop.txt) \
    | cut -f 1,2,3,7,11,12 | gzip > transaction.txt.gz

こうやって商品マスター、店舗マスター、トランザクションの三つに分割しています。このcsvはたまたま要素内に半角カンマがないことがわかっていたので、切り出しは単純にカンマ区切りでOKでしたが、ダブルクオーテーションの中にカンマや改行が入っていたらもっと厄介なことになっていました。カラム3と7がそれぞれ商品マスタと店舗マスタの固有IDです。pos.csv.gzはめちゃめちゃ巨大なので、展開したデータを置いておきたくないですし、さらに文字コードを変換したデータを置くのも厳しいので、展開と文字コード変換を一度だけ実施した上で三つの作業を同時にやりたかった、というわけです。

最近、マーケティング系の研究会でヒットとクチコミとの関係モデルを研究していまして（参考：電通/吉田財団助成研究から、画期的な成果（「環メディア」の発見と「コクーンブレークモデル」の実証）が生まれた（法政大小川教授のブログ））、日頃からヒット商品とそのクチコミは興味深くウォッチしているのですが、今年は映画が非常に面白いですね。「聲の形」、「シン・ゴジラ」、「君の名は。」、そして「この世界の片隅に」。どれもネットメディアやSNSでかなり話題になっており、興行成績もそれに合わせてこれまでの映画と異なる推移を示しているようです。

その中でも大注目しているのが、「この世界の片隅に」です。

上のリンクで示したモデルは、草の根から発信が始まって、ある時期に中規模ネットメディアに見つかり、それが拡散を広範化させ、それを繰り返すうちにテレビに見つかって大ブレークするというメカニズムを観察しています。最初のステップは共通の趣味の人が訪れるまとめサイト、次は中規模のネットニュース、それがYahoo!JAPANのニュースに繋がり、地上波ワイドショーで取り上げられてブレークするのが共通したヒットの法則でした。しかしこれを裏返すと、大ヒットの裏には必ずテレビの影響があるということ。大ヒットにはテレビは欠かせないメディアであり、テレビの影響なくして大ヒットには繋がらない（もしくは、非常に難しい）というのが、コクーンブレークモデルの結論の一つです。

しかし、映画「この世界の片隅に」は、もしかしたらそのモデルを覆すのかもしれません。なぜなら、これまでテレビでのプロモーションがほとんど行われていないからです。さらには、SNSやネットメディアでこれほどまでに盛り上がっているのにもかかわらず、ワイドショー等が取り上げない（上のコクーンブレークモデルでは、ある程度の盛り上がりがあったら、テレビがそれを見つけて再拡散するというメカニズムが働く）という希有な事例でもあります。

プロモーションが行われない理由、また、テレビがニュースとして取り上げない理由はここでは関知するところではありませんし、議論する気もありません。私が科学として興味があるのは、テレビ露出がこれまでほとんどないという事実に対して、そのような商材がなぜこれほどまでに情報拡散され、消費されているのか、その情報がどのような拡散のメカニズムを経由したのかというところです。特に、映画館ではお年寄りの姿も多いとのことですが、しかし、お年寄りはネットを利用するよりも、テレビから情報を得るのが圧倒的というのがこれまでのモデルでした。もしかしたら、いわゆるM3F3層のネット利用率、ネット使いこなし率は、現在思っているよりもかなり高いのかもしれない、などの興味がどんどん湧いて出てきます。もちろん、実際のデータを見てみないことには本当のところはわかりませんが、想像するだけでも、これまでのプロモーションの通説を覆すいろんな要素が、この映画の拡散プロセスには詰まっているような気がします。


ところで。

本作品、もちろん科学としての興味も大きいのですが、私はこの映画自体も素晴らしいと感じましたし、作品の素晴らしさが突き抜けているからこそ、大規模プロモーション無しでこれほどの拡散が観察されるのだろうと思いました。実は二度ほど観に行きまして。私の祖母がちょうどすずさんと同じくらいの年齢だと思います。幼い父を抱えて空襲の中を逃げたという当時の話を聞いていましたから、あの光景が非常に身近に感じられます。防空壕の中にいて、落ちてくる爆弾の音のリアリティも鳥肌がたつほどの恐怖をかき立てます。そして、クライマックスからエンディングにかけての子持ちの親としての感慨もあります。つまり、そもそも作品として素晴らしいし、また、観終わった後、これを身近な人にもぜひ見てもらいたいと思う気持ちが非常に強く感じられるのが、本作品の特徴なのだろうと思います。

最初はわずか63館で始まったとのことですが、年明けには190もの映画館での上映が決まっているとのこと。また、海外上映への動きも始まっているとのことで、今後も興味深く観察していきたいと思います。あとは、データがあったらぜひ分析してみたいですね。

早いもので、新会社を設立してから三ヶ月が経ちました。ブログの引越がずるずると先延ばしになっていたのにやっと手を付けまして、全部こちらの jinya.seesaa.net に引っ越しました。

設立したのは「株式会社ゴーガ解析コンサルティング」で、以前に取締役を務めていました株式会社ゴーガからデータ分析・解析のチームを子会社で分社化しました。旧社で10年やったので、もう一度経営方針を考え直して心機一転といったところです。目論見は二つありまして、一つは解析だけのチームで子会社化することで、データ分析コンサルティングの分野でより一層の機動力を発揮することなのですが、もう一つは人事面の見直しです。

ウェブの世界とデータ分析の世界は似ているようでちょっと違っていて、ウェブの世界は習うより慣れろでぐいぐい開発の経験値を溜めて一人前になる印象がありますけれど、データ分析の世界は地味な作業や勉強とコミュニケーションを積み上げなければなかなか一人前にはなれません。技術の一点突破ではなくて、技術と理論とコミュニケーションを偏りなく高めていく必要があります。少なくとも、私が志向するデータ解析コンサルティングとはそういう方向性です。その結果、どうしても人材育成には時間がかかる。その時間差みたいなものを解消しようというのが、今回の分社化のもう一つの目論見です。

数年前から「ビッグデータ」と共に「データサイエンティスト」なるキーワードがもてはやされて、業界団体ができたり（弊社も参加しています）人材争奪合戦が激しさを増してきたりしていますが、サイエンスは流行り廃りに左右されるものではありません。数学や統計学などの知識と経験をしっかりと基礎から積み上げて、そこに最新のデータエンジニアリングのテクニックと、ビジネスの知見を織り交ぜてコンサルティングができる、そういう人材を育成していくことが、私の次の10年の取り組みです。

最近、自分の計算用マシンにWindows10のアニバーサリーアップグレードを入れてから、調子が悪いなぁと思っていたのですが、メモリの使用可能域が半分になっていることに気がつきました。[システム]を見ると、実装メモリ(RAM)は16GBなのに、その横に括弧書きで（使用可能 7.68GB）と書いてあることを発見。リソースモニターなどもよくよく見ると最大が8GBになっていて、ハードウェア予約済みに半分の8GBが持って行かれている状況。

ググると、BIOSだ、内蔵グラフィックだ、ブートだ、といろんな説が出てきたのですが、どれを実行しても効果がなかったり、そもそも項目がなかったり。と、その中の一つの解決方法に、「メモリを刺し直す」っていうのがあって、まさかそんなとは思いつつ最後の手段として実施してみたところ、見事解決。そんなバカな、と思いながらも、やってみるものですね。

ところで、アニバーサリーアップグレードはもちろん、Ubuntu on Windowsを試すために入れたのですけれど、これは開発者モードをONにすると自動的にsshもONになります。それなりに危ないので、UoWを使う際はsshのサービスは無効化しておくのがよいと思います。

「株式会社ゴーガ解析コンサルティング」を設立します。

ニュースリリース：
http://www.goga.co.jp/company/news20160523.html

営業開始は7月1日より。ただいま鋭意準備中です。
今後ともよろしくお願いします。

MS Office 2016の配布が始まりましたので、さっそく更新して使ってみたところ、ワードで見慣れない編集記号が。Pの鏡文字みたいな圧迫感のある記号がすべての改行（もしくは改段落）のところに出現しました。

編集記号は表示させたいけれど、見た目がこれでは気持ちよく編集することができない、どうしたものか・・・と思って解決策をググっていたところ、さすがネットの世界、既に解決なさっている方がいらっしゃいましたので、そちらをご案内します。

→ 「Word 2016 for Windowsの段落記号が変だった話（らむれーずんめも）」（リンク先修正しました）

問題は「規定の言語」が日本語から英語に変わってしまっていたことにあったようです。MSではよくあることですね（笑）。

RAID1を使ってHDDの容量を増やした話です。OSはWindows10, マシンは普通のデスクトップ（EPSONのMR4000）ですが、「RAIDキット」を搭載しています。おそらくRAIDカードを挿しているのだろうと思いますが、詳しくは知りません。

RAIDはRAID1を利用していて、片方のHDDが壊れたので換装しました。本当ならば同じメーカー、同じ型番のHDDを購入して片方だけを取り替えればいいのですが、ちょっと古いマシンだったので同じモノは見つからず。（実際には、中古または新品でもマニア向けな価格ならばありましたが、500GBのHDD一本に4万円は出せない！）そこで、普通に販売しているHDDを二本購入し、両方を入れ替えることを試みました。

破損した一本をまず取り替え、リビルドした後、まだ生きている方も交換してリビルド。特に問題も発生せず、シンプルに二本のHDDの換装完了。

ところで、古いHDDは500GBで、新しく購入したHDDは1TB。そんなに大量にディスクを使う方ではないので500GBでも十分なのですけれど、1TBが安かったので、残りの部分は捨てるつもりで1TBのHDDを選びました。しかし実際に搭載してみるとやはり500GBしか認識されないのはどこか勿体ない…と思っていたところ、RAID容量を増加させるメニューをRST（Intel Rapid Storage Technology）の中に発見しました。

昔は、RAID1である容量のボリュームを一度構成したら、あとはその容量にずっと縛られるのが普通だったと思います。つまり、最初に500GBで構成すると、あとはいくら大きい容量のディスクに入れ替えようとも、そのボリュームの容量は500GBから変えられなかったと思います。実際、適当にググると、容量を増やしたい！という記事は見つかるのですが、できませんという返事ばかりが引っかかりますね。

ところが時代は変わるもので、上述の Intel Rapid Storage Technology にはRAID1ボリュームの容量を増やす仕組みが搭載されていました。RSTのボリュームの管理の画面に行くと、ボリュームよりも両方のHDD容量が大きい場合にそのリンクが出現しますので、そこから容量を拡大できます。但し、初期化にかなり時間がかかります。

RSTで初期化が終われば、再起動の後、増えた分がディスクの未割り当て領域として[ディスクの管理]などに出現します。あとはその領域を別パーティションとしてマウントするなり、既存のパーティションを拡張するなりすればよいです。

ちなみに、私のケースではCドライブがベーシックディスクでしたので、[ディスクの管理]ではCドライブを拡張することができませんでした。ですが、diskpartならばできるということで、コマンドプロンプトからdiskpartでdiskpartを起動、selsct disk, select partition, select volumeでCドライブを選び、extendで拡張し、これでCドライブが500GBから1TBになりました。



なお、最後に定型文ですが、この通りやっても容量を増やせる保証はありませんので、お試しになる方は自己責任で。

自宅での話。

次男坊が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が出てきました。

回帰モデルのギブスサンプラーを導出します。教科書は「マルコフ連鎖モンテカルロ法」（豊田秀樹著）。こちら 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\)を必要としないこともわかる。

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\)に依存しないと言うことは、単純にギブスサンプリングでの計算量が半分ほどに減るわけで、そもそも計算量を減らすためにメトロポリス-ヘイスティングスではなくギブスサンプラーを使っていることを考えると、効率がより向上することになります。

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