EbiMaru’s diary

ラマン分光とか機械学習とか研究している京大生のブログ

インターンで論文を出すことになった話

 M1の夏休みに某大企業の研究開発を担当している所(以下ONSHA)に3週間インターンに行った感想とかを書きます.

インターン前まで

 そもそもなんでONSHAに限らずインターンシップに参加しようとしたかというと,一応業界研究はしておきたかったし,まあ就活にも悪い影響は及ぼさないでしょと思ったからです.特に深い意味は考えていなかったですね.いやこのひと,普通にマ○ナビとかリ○ナビにいいように利用される資本家の犬でしかないやん...ちなみにですが,選考は1勝6敗でしたので,参加するとしたらONSHAしかありませんでした().
 正直ONSHAは一番選考が厳しそうだと思っていたので,通ったときから北斗の拳的世紀末を想像していました.というのもインターンは即戦力になるようなスキルをもっていないと通らないと思っており,ONSAHのインターンのテーマを見た時,自分の研究にもスキルにもドンピシャで合致するテーマはなく,「あ,ここには通らないな」と思いました.そこで,第1志望はせめて戦えそうな機械材料系のものにして,第2志望以下は機械学習系のものを選択しました.しかし,来たるべき7月,選考の結果は「第2志望:機械学習ヲ用ヒでーたノ解析ヲセヨ」でした.震えました.周りはどんな強者がいるのだろうと.自分は果たして,スキルをチョット盛った罪悪感と劣等感に耐えられるのかと.でも私の中のクズ要素は,この不安を見事に押しのけ,インターンの1週間前まで開発に使うと言われていたPyTorchを触らせず,PRMLも読めと言われていた論文も読ませなかったのです...!まあ最終的には簡単な回帰の実装をして,論文も読みましたが...あとパワハラにあった場合に備えて,レコードできるペン買いました.要りませんでしたが.

業務内容編

 私に与えられ課題は,まあ画像処理の一種で,業務自体はある意味で単調でした.論文を読む→git clone and 実装→モデルを学習→結果の吟味→議論の繰り返しです.僕が配属された部署は,「機械学習を用いて現場の課題解決」とか「機械学習で金儲け」ではなく専ら論文のための研究を行っていたようです.途中で気づいていましたが,「これは大学と大して変わらんな」と.企業研究という点ではONSHAはあまり適切ではなかったかもしれません.他のインターン生の話を聞くと確かに大学の研究と似てはいるが,大学よりも期日などをしっかり定め限りある時間を最大限に利用とすることを教えられる人が多かったらしく,一方僕は特に一日のノルマなどの定義もなかったのでのびのびと研究していました.まあ,業務中はほとんどTwitterできませんでしたが().同じチームの方々は,皆わからないことがあればフランクに教えてくれ,間違ったことをいえば特に嫌味な風でもなく指摘してくれました.
 さて,自分はどのように活躍したかというと,そんなには手出しできなかったですね.一応自分からアイデアをポンポン提案とかできればよかったのですが,言われた手法を実装することしかできなくてちょっと悔しかったですね.僕は実はNNのモジュールを触るのは今回が初めてでした.ML系の論文も初めて読みました.まあ,そんな素人同然...!な状況で,勝てるのは赤木しげるくらいなわけで.
 そういった初めての経験をしていくなかで,思ったのは少なくとも機械学習を生業とし,NN職人として生きることだけはやめようということです.レッドオーシャンすぎる.その上これは泥沼だと.落合陽一氏が言ってたのですが,「デジタルネイチャーの学生にわかって欲しいことはみんなが憧れていて練習したり努力すればやがてなれそうな気がする職業っていうのはレッドオーシャンだってことだよ.」というのが刺さります.まあ,もとよりその気もなかったのですが.しかし機械学習は今後も発展し,データ解析の点で新たな知見を私達に運んでくれるのだと信じています.だから,最先端を追うことはやめてはならないとも思いました.私達実験屋はフロンティアを開発してくれた勇気ある人たちの開拓してくれた土地を,ありがたく豊かにすることを考えるべきだと,改めて思いました.

社風編

 まあ,予想はしていましたが,職場はホワイトでした.残業はほぼ0.ある人は残業時間がマイナスになったと言っていました.育休・有給の類は簡単に取れ,「気分が乗らなかったら年休とる」とかも普通だったようです.いや大学かよ.それだけでなく,各人が大人の都合をわきまえ他人を最低限尊重し,無意味に他人を蹴落とすようなことのない職場でした.まあ,職場のシステム的に出世や権力といった概念がほとんど役に立たないところだったので当たり前といえば当たり前なのかもしれません.少し話は逸れますが,残業は悪ですね.労働者が好きでやっているならわかりますが,残業込みで働かないといけないのは監督者の無能と怠惰ですよ.21時まで残って家まで着いて色々済ませ,朝には眠気眼に発破をかけるとか,自分にできたものじゃないと思います.
 職場の人達はさっきも申したとおり,いい人たちでした.その上ドチャクソ頭がいいのです...!機械学習の研究をしていたわけですが,皆,機械学習系のトップカンファレンスに論文を通すことを目標に研究していたようです.物理・情報・統計・材料・化学・生物.これらをまんべんなく,しかし各々の知識が連携できるくらいには深く理解していたように感じました.よく「人はスペシャリストとジェネラリストどちらを目指すべきか」みたいな議論がありますが,その点で言えばONSHAの方々はジェネラリストに近かったと思います.でも,ある程度論文が書けたり,設計開発ができたりしなければジェネラルでも役にたたないのだと気付かされました.なんとしても,そこまでたどり着かなければ.
 ただ,ONSHAの社員の方に言われたことで「ここは大学と大して変わらない.研究はできても,開発には遠い分野の人が多い.研究ができるだけでは直接は金は生まれない.具体的には例えばうちに入った後に転職からの中途採用では不利になる可能性がある.しかも取り潰しかどうかはお上の判断が全て.経営が傾いたら真っ先に切られる存在.実際企業の基礎研究部門はバブルの頃に流行ってまだうちは生き残っているだけ.いわばバブルの化石」というのもありました.研究は道楽なのか.ONSHAレベルの所でもそのような声が聞こえるのは悲しいとは思いますが,それは一つの現実なんでしょう.企業人として生きるには,やはりマネジメント・ビジネス感覚は必須なのだと思い知らされました.好きなことで食っていくのは難しいが,得意なことで食っていくのは簡単だというのは真に感じます.
 それでも,ONSHAに惹かれるのは単純に職場の雰囲気とかだけではなく,人間的にな魅力の多い人がいるからです.まあサンプルはONSHAしかありませんが.ある人は,ONSHAの仕事以外に,ボランティアで障害者や外国の貧困層の支援活動を行っていました.自分が他人から攻撃されない領域で他人を叩いてヘイトを撒き散らしている間に,真に弱者の立場に立とうとして現に支援をしているひとがいる.この事実は辛いっすねー.自分が情けなく思えちゃいますよー.過去に拘泥するはいい加減やめたい.実際どうしようもないこともありますが.そうでない人も,大体Ph.D持ちで,なんというか人生を自分より10倍くらい経験してる感がありましたね.あ,あと30過ぎてる人はほとんど既婚者子持ちでした.

研究の結果

 まあまあうまく行きましたよ.オープンなデータセットを解析していて,もともとは別のタスクを行っていたのですが,分類問題にも応用したら元論文の精度89%だったのが,96%に向上しました.ワールドレコードが97%とかだったので,個人的にはかなり意義の高い研究だったのではないかと考えています.この成果でイコール・ コントリビューションで論文を出す話がまとまってきています.まあ,僕は言われたこと実装しただけですが(爆).最終的にどうなるかは神のみぞ知るです().

最後に

 実際の生産業務に関わりをもつひとたちに比べれば多少給料は安いものの,私のいたところは上記のようなホワイト度で,周りの人間も尊敬できる人が多かったと思うので,今回のインターンで行きたい会社の上位に食い込みましたね.しかし,ONSHAに限らず,閉鎖的な空間にのみ居ては,人間的にいずれ死んでいってしまうのではないかという危機感は覚えました.そのためだけに転職などをするというのはドMの所業だとは思いますが,これからの時代に評価されるキャリアとは何かを,学生のうちに考えておきたいものです.あとビジネス感覚身につけたいな.とにかく自分の見ている世界は狭いなと思いました.たくさん国際学会に行って,いろいろな人とおしゃべりしないとですね.
 ONSHAのインターンシップスキルアップにはなったし,すごい人もいるのだと感じれたから,その点は☆5つですね.ただ,就活自体にはそんなに経歴にプラスになるようなことはなかったと思うので,そこはもう腹をくくるしかない.人生経験としてはとても良かったですかね.アカデミックとしての成果も得られそうだし.まあ就活の参考にはあまりならない駄文だったかもしれませんが,少しでも情報を提供できたらいいですわ.

f:id:EbiMaru:20190831224600j:plain
インターン中の休日に訪れた伊勢神宮の内宮の入り口

f:id:EbiMaru:20190831224511j:plain
喫茶コンパルで初めて食べたモーニング

調和振動子の古典・量子

そうだ.量子論で一次元調和振動子を解こう

 というわけで確率密度関数を求めます.アニメはpyてょnしか書けないので,それで作ります.

古典力学

調和振動子ハミルトニアンは近所の佐藤少年でも知ってる

\displaystyle{
H = \frac{p^2}{2m} + \frac{m \omega^2}{2}x^2.
}

(m:質量,\omega:振動数).解き方は幼稚園児でもわかるので省略.

\displaystyle{
x(t) = x(0) \cos \omega t + \frac{p(0)}{m \omega} \sin \omega t
}
\displaystyle{
p(t) = -m \omega x(0) \sin \omega t + p(0)\cos \omega t
}

量子力学

面倒くさいなあ.まずは状態ケットの時間変化を見るか(Schrödinger描像にしときますわ.あと表記の都合上,波動関数が時間に依存するかどうかは引数の中身で区別しています).

\displaystyle{
H \left| n \right> = \left( \frac{1}{2} +n  \right) \hbar \omega \left| n \right>
}

だから,

\displaystyle{
\mathcal{U(t)}  \left| n \right> =  \exp \left[ -i  \frac{H t}{\hbar}\right] \left| n \right>= \exp \left[ -i  \left( \frac{1}{2} + n \right) \omega t \right] \left| n \right>
}

(\mathcal{U}(t):時間発展演算子)です.エネルギー固有関数が求まれば,波動関数も求まることがわかります. \sigma = \sqrt{\hbar/ m \omega}を長さの規格化定数とします.それらを求める前にどのくらいの長さか考えてみましょう.水素分子のばね定数が510 N/m程度

http://www.comp.tmu.ac.jp/qc/toudai/2007-11-9.pdf, p. 6.

mは換算質量(m = m_\mathrm{hydrogen}/2)であることを考慮して水素分子の振動数は

\displaystyle{
\omega \approx \sqrt{510/(1.674 \times 10^{-27} / 2)} = 7.81 \times 10^{14} \hspace{5pt} \mathrm{s}^{-1}
}

程度です.また,\sigma

\displaystyle{
\sigma = \sqrt{\hbar/\sqrt{km}}  \approx  \sqrt{1.055 \times 10^{-34} / \sqrt{1.674\times10^{-27} \times 510 /2}} 
}
\displaystyle{
= 1.27 \times 10^{-11} \hspace{5pt} \mathrm{m} 
}

で,大体0.1 Åだから,そんなものかあといった感じです.
 エネルギー固有関数と波動関数を求めていきましょう.u_n(x')を位置表示のエネルギー固有関数,v_n(p')を運動量表示のエネルギー固有関数とします.また,\psi_n(x',t)を位置表示の波動関数\phi_n(p',t)を運動量表示の波動関数とします.

\displaystyle{
\sqrt{\frac{m \omega}{2 \hbar}} \left< x' \left| a \right|0 \right> = \sqrt{\frac{m \omega}{2 \hbar}} \left< x' \left| \left( x + \frac{ip}{m \omega} \right) \right|0 \right> = 0
}

( a は消滅演算子)から

\displaystyle{
\left( x' + \sigma^2 \frac{d}{dx'} \right) \left< x' | 0 \right> = 0
}

この微分方程式を解くと基底状態のエネルギー固有関数が求められ

\displaystyle{
u_0(x') = \left< x' | 0 \right> = \left( \frac{1}{\pi \sigma^2 } \right)^{1/4} \exp \left[ - \frac{1}{2} \left( \frac{x'}{\sigma}\right)^2 \right]
}

です. また

\displaystyle{
   \left< x' | 1 \right> = \left< x' | a^\dagger| 0 \right>
}

から

\displaystyle{
  u_1(x') = \left< x' | 1 \right> =\left( \frac{1}{\sqrt{2} \sigma } \right) \left( x' - \sigma^2 \frac{d}{dx'} \right) \left< x' | 0 \right> = \left( \frac{4}{\pi \sigma^2} \right)^{1/4} \frac{x'}{\sigma} \exp \left[ - \frac{1}{2} \left( \frac{x'}{\sigma} \right)^2 \right]
}

まあ,エルミート多項式使ったほうが早いんですけどね.さて,波動関数がもとまったので粒子の存在確率\rho_xを求めます.初期状態が重ね合わせ状態じゃないと振動しないので

\displaystyle{
   \left| \alpha;t=0 \right> = \frac{1}{\sqrt{2}} \left| 0 \right> + \frac{1}{\sqrt{2}} \left| 1 \right>
}

とします.ゴリゴリ計算して,

\displaystyle{
   \rho_x(x',t) = \left| \frac{1}{\sqrt{2}} \psi_0(x',t) + \frac{1}{\sqrt{2}} \psi_1 (x',t) \right|^2 
}
\displaystyle{
=  \frac{1}{ \sigma \sqrt{\pi}} \ \exp \left[ - \left( \frac{x'}{\sigma} \right)^2 \right] \left[  \frac{1}{2}  +  \left( \frac{x'}{\sigma} \right)^2 + \sqrt{2} \frac{x'}{\sigma} \cos \omega t \right]
}

です.確かにx'積分すれば,1になりますね.次は運動量の確率密度\rho_p(p,t)を求めましょう.運動量表示のエネルギー固有関数v_0(p')v_1(p')を求めます.簡単に計算できて(←言ってみたかっただけ)

\displaystyle{
v_0 (p') = \frac{1}{\sqrt{2 \pi \hbar}} \int dx' \exp \left( \frac{- ip' x'}{\hbar} \right) \psi_0 (x')
}
\displaystyle{
= \sqrt{ \frac{\sigma}{ \hbar \sqrt{\pi}}} \exp \left[ - \frac{1}{2} \left( \frac{p' \sigma}{\hbar} \right)^2 \right]
}

v_1(p')は少し面倒ですが,

\displaystyle{
v_1 (p') = \frac{1}{\sqrt{2 \pi \hbar}} \int dx' \exp \left( \frac{- ip' x'}{\hbar} \right) \psi_1 (x')
}
\displaystyle{
= \frac{1}{\sqrt{2 \pi \hbar}} \int dx' \left[  \exp \left\{ - \frac{1}{2} \left( \frac{x'}{\sigma} + \frac{i p' \sigma}{\hbar} \right)^2 - \frac{1}{2} \left( \frac{p' \sigma}{\hbar} \right)^2 \right\}  \left( \frac{x'}{\sigma} + \frac{i p' \sigma}{\hbar} \right) - \frac{ip'\sigma}{\hbar} \exp \left\{ - \frac{1}{2} \left( \frac{x'}{\sigma} + \frac{i p' \sigma}{\hbar} \right)^2 - \frac{1}{2} \left( \frac{p' \sigma}{\hbar} \right)^2 \right\} \right]
}
\displaystyle{
= \sqrt{ \frac{2 \sigma}{ \hbar \sqrt{\pi}}} \frac{- i p' \sigma}{\hbar}\exp \left[ - \frac{1}{2} \left( \frac{p' \sigma}{\hbar} \right)^2 \right]
}

なので,運動量表示の波動関数も求められます.確率密度関数\rho_p(p',t)

\displaystyle{
   \rho_p(p',t) = \left| \frac{1}{\sqrt{2}} \phi_0(p',t) + \frac{1}{\sqrt{2}} \phi_1 (p',t) \right|^2 
}
\displaystyle{
  = \frac{\sigma}{\hbar \sqrt{\pi}} \exp \left[ - \left( \frac{p' \sigma}{\hbar} \right)^2 \right] \left[ \frac{1}{2} + \left( \frac{p'\sigma}{\hbar} \right)^2 - \sqrt{2} \frac{p'\sigma}{\hbar}  \sin \omega t  \right] 
}

割と対称的できれいですわね(もっと簡単に求める方法があったら教えてください><) さて,アニメーションだ. \hbar=\sigma=1,\omega = 2 \piとします.結果は以下のとおりです
f:id:EbiMaru:20190721204339g:plain
ミュー空間(と言っていいのかな...)でのカラーマップはこんな感じ
f:id:EbiMaru:20190721204351g:plain
なんかあってる気がしないな笑.期待値はちゃんと楕円軌道になっていると信じたい.
 この調子で高エネルギーの初期条件の運動も求めていきます.さて,ここであらかじめ作っておいた準位0~3の固有関数を用意します(笑). 計算の都合上,X' = x/\sigmaP' = p \sigma /\hbar として,

\displaystyle{
u_0(X') = \left( \frac{1}{\pi \sigma^2 } \right)^{1/4} \exp \left[ - \frac{1}{2} X'^2 \right]
}
\displaystyle{
  u_1(X') = \left( \frac{4}{\pi \sigma^2} \right)^{1/4} X' \exp \left[ - \frac{1}{2} X'^2 \right]
}
\displaystyle{
  u_2(X') = \left( \frac{4}{\pi \sigma^2} \right)^{1/4} \left[ X'^2 - \frac{1}{2} \right] \exp \left[ - \frac{1}{2} X'^2 \right]
}
\displaystyle{
  u_3(X') = \left( \frac{4}{9 \pi \sigma^2} \right)^{1/4} \left[ X'^3 - \frac{3}{2} X' \right] \exp \left[ - \frac{1}{2} X'^2 \right]
}

積分

\displaystyle{
\int dX' \exp \left[ - \frac{1}{2} \{ X' + iP' \}^2 - \frac{1}{2} P'^2 \right] = \sqrt{2 \pi} \exp \left[ - \frac{1}{2} P'^2\right] 
}
\displaystyle{
\int dX' \exp \left[ - \frac{1}{2} \{ X' + iP' \}^2 - \frac{1}{2} P'^2 \right] X' = i \sqrt{2 \pi}  \exp \left[ - \frac{1}{2} P'^2\right] \left( - P' \right)
}
\displaystyle{
\int dX' \exp \left[ - \frac{1}{2} \{ X' + iP' \}^2 - \frac{1}{2} P'^2 \right] X'^2 = \sqrt{2 \pi}  \exp \left[ - \frac{1}{2} P'^2\right] \left( - P'^2 + 1 \right)
}
\displaystyle{
\int dX' \exp \left[ - \frac{1}{2} \{ X' + iP' \}^2 - \frac{1}{2} P'^2 \right] X'^3 =i \sqrt{2 \pi}  \exp \left[ - \frac{1}{2} P'^2\right] \left(  P'^3 - 3P'  \right)
}

を用いてこれに対する運動量表示の固有関数を求め

\displaystyle{
v_0 (P') = \sqrt{ \frac{\sigma}{ \hbar \sqrt{\pi}}} \exp \left[ - \frac{1}{2} P'^2 \right]
}
\displaystyle{
v_1(P') = i \sqrt{ \frac{2 \sigma}{ \hbar \sqrt{\pi}}} (-P')\exp \left[ - \frac{1}{2} P'^2 \right]
}
\displaystyle{
v_2(P') = \sqrt{ \frac{2 \sigma}{ \hbar \sqrt{\pi}}} \left( - P'^2 + \frac{1}{2}  \right) \exp \left[ - \frac{1}{2} P'^2 \right]
}
\displaystyle{
v_3(P') = i \sqrt{ \frac{2 \sigma}{ 3 \hbar \sqrt{\pi}}} \left(  P'^3  -\frac{3}{2} P' \right) \exp \left[ - \frac{1}{2} P'^2 \right]
}

さて,同じ重みで重ね合わせされた,準位が1つだけ違う状態の重ね合わせの初期状態での確率密度関数

\displaystyle{
\rho_{X}^{(12)} (X',t) =  \frac{2}{ \sigma \sqrt{\pi}} \ \exp \left(- X'^2 \right) \left[  X'^2 +  \left(  X'^2 - \frac{1}{2} \right)^2 + 2X' \left(  X'^2 - \frac{1}{2} \right) \cos \omega t \right]
}
\displaystyle{
\rho_{X}^{(23)} (X',t) =  \frac{2}{ \sigma \sqrt{\pi}} \ \exp \left(- X'^2 \right) \left[ \left(  X'^2 - \frac{1}{2} \right)^2 + \frac{1}{3} X'^2 \left( X'^2 -\frac{3}{2} \right)^2 +  \frac{2}{\sqrt{3}}X' \left(  X'^2 - \frac{1}{2} \right) \left( X'^2 - \frac{3}{2} \right) \cos \omega t \right]
}
\displaystyle{
\rho_{P}^{(12)} (P',t) =  \frac{2 \sigma}{\hbar \sqrt{\pi}} \exp \left( - P'^2 \right) \left[ P'^2 + \left( P'^2 - \frac{1}{2}\right)^2 - P' \left( P'^2 - \frac{1}{2} \right)^2 \sin \omega t \right]
}
\displaystyle{
\rho_{P}^{(23)} (P',t) =  \frac{2 \sigma}{\hbar \sqrt{\pi}} \exp \left( - P'^2 \right) \left[ \left(  P'^2 - \frac{1}{2} \right)^2 + \frac{1}{3} P'^2 \left( P'^2 -\frac{3}{2} \right)^2 -  \frac{2}{\sqrt{3}}P' \left(  P'^2 - \frac{1}{2} \right) \left( P'^2 - \frac{3}{2} \right) \sin \omega t \right]
}

アニメにするとこんな感じです.

f:id:EbiMaru:20190721204410g:plain
確率密度関数(初期状態:1 and 2)

f:id:EbiMaru:20190721204428g:plain
確率密度関数のカラーマップ(初期状態:1 and 2)

f:id:EbiMaru:20190721204446g:plain
確率密度関数(初期状態:2 and 3)

f:id:EbiMaru:20190721204457g:plain
確率密度関数のカラーマップ(初期状態:2 and 3)

やはり初期エネルギーが大きいと,節が増えていくんですね.
 次に,期待値と偏差を考えていきます.

\displaystyle{
 a = \frac{1}{\sqrt{2}} \left( X + iP \right), \quad a^\dagger = \frac{1}{\sqrt{2}} \left( X - iP \right)
}

で,

\displaystyle{
 X = \frac{1}{\sqrt{2}} \left(  a^\dagger + a  \right),\quad  P = \frac{i}{\sqrt{2}} \left( a^\dagger - a \right)
}

です,初期状態を\left| \alpha;t=0 \right> = \left( \left| n \right> + \left| n+1 \right> \right)/\sqrt{2}とすると

\displaystyle{
 \left| \alpha;t \right> = \frac{1}{\sqrt{2}}\left( \exp \left[ - \frac{1}{2} (n + 1) \omega t \right] \left| n \right> + \exp \left[ - \frac{1}{2} (n + 2) \omega t \right] \left| n+1 \right> \right)
}

から,

\displaystyle{
 \left< X \right> = \sqrt{\frac{n+1}{2}} \cos \omega t , \quad \left< P \right> = -\sqrt{\frac{n+1}{2}} \sin \omega t
}

また,

\displaystyle{
 X^2 = \frac{1}{2} \left( a + {a^\dagger}^2 + a a^\dagger + a^\dagger a \right), \quad P^2 = \frac{1}{2} \left( a a^\dagger + a^\dagger a - a^2 + {a^\dagger}^2 \right)
}

なので,

\displaystyle{
 \left< X^2 \right>= n + 1 ,  \quad \left< P^2 \right> =  n + 1
}

となり,偏差は

\displaystyle{
 \Delta X = \frac{1}{2}\sqrt{(n+1) (3 - \cos 2 \omega t )}, \quad \Delta P =  \frac{1}{2}\sqrt{(n+1) (3 + \cos 2 \omega t )}
}

です.xの広がりは\sigma \propto 1/\sqrt{mk} に,pの広がりは\sigma^{-1}\propto\sqrt{mk}に比例し,いずれも\hbarに比例するので,古典極限は案の定,\hbar \rightarrow 0であり,ばね定数に対する質量の比が小さいほど位置が,大きいほど運動量が不確かになることが分かります.

具体例で学ぶベイズモデル選択

ベイズ的推定におけるモデルの選択方法を書いていきます.基本的にはベイズ的線形回帰を例に取り,基底の数や種類を変えて,どのような基底をどのくらい用意するべきかを書いていきます.ここからエビデンス近似なんかに繋がって行きます.

Keywords

ベイズ統計,ベイズモデリング,モデル選択

 モデルとは?

モデルというだけだと抽象的ですよね.例えば,事前分布であったり,ハイパーパラメータの値だったり,仮定する確率分布だったり,基底関数だったり...これらは全部モデルなんですが,それらをまとめてビショップ先生に敬意を表して(?)\mathcal{M}_i と書いておきます(ただし今回いじるのは基底の数と種類だけです).ハイパーパラメータなんかは交差検証でそれなりの値を出せたりします(こちれはどちらかというと頻度分布主義的な手法で今回は使いません).更に,赤池情報基準(AIC),ベイズ情報基準(BIC)とかWAICなんてものもありますね.これらに共通するモデル選択の目標とは,得られたデータに過度にフィットしないようにしつつ(つまり過学習を避ける),単純すぎないで予測精度が悪くならないようなモデルを選ぶ(つまり観測されたモデルパラメータ\check{\boldsymbol{w}}と推定された\hat{\boldsymbol{w}}の差を小さくする)ことです.今回はハイパーパラメータの大小の比較ではなく,xの累乗と別の「表現力のいい」関数を用いて比較していきます.ハイパーパラメータに関しても同様の思想で「エビデンス近似」と呼ばれる手法で最適化できます.
 \sin 2\pi x+1にガウシアンノイズが乗っているような系で考えていきます.基底としてx^{m}(m=0,1,..,M-1)を使うとして,何個使えばいいでしょうか?あるいは,得られたデータの分布から\sin 2 \pi x1で近似したいと思ったときどちらがいいでしょうか?現実的には,真の分布なんてわからないわけですから,「最もそれっぽい」ものを選ぶ必要があります.前準備としてベイズ的線形回帰の式を結果だけ示します.

ベイズ的線形回帰

ベイズ的線形回帰とは,実現値の集合\mathcal{D} = \{ \check{\boldsymbol{t}},\check{\boldsymbol{x}} \}が与えられたとき,再構成の目標値tを,M個の基底関数\phi_m(x)を用い,それにかかる係数\boldsymbol{w}=(w_0, w_2,...,w_{M-1})^\mathrm{T}から推定する手法で,例えば事前分布を平均\boldsymbol{m}_0,共分散を(\alpha \beta)^{-1} \boldsymbol{S}_0の多変量正規分布と仮定し,尤度を分散\sqrt{1/\beta}の多変量正規分布とすると,係数に関する事後分布は正規分布となり


p(\boldsymbol{w}|\mathcal{D}) = \mathcal{N} (\boldsymbol{w}|\boldsymbol{m}, \beta^{-1} \boldsymbol{S})\\
\boldsymbol{S}^{-1} = \alpha \boldsymbol{S}_0^{-1} +  \boldsymbol{\Phi}^\mathrm{T} \boldsymbol{\Phi} \\
\boldsymbol{m} = \boldsymbol{S}(\boldsymbol{\Phi}\check{\boldsymbol{t}}+\alpha \boldsymbol{S}_0^{-1} \boldsymbol{m}_0)

です.ただし\boldsymbol{\Phi}は計画行列(ref. PRML p.139)です(事後分布の平均が\alphaのみに依存するように事前分布の分散を調節しています.\alphaは事前分布の影響の強さを示します).確率モデルp(t|x,\boldsymbol{w},\beta)

\displaystyle{
p(t|x,\boldsymbol{w},\beta) = \mathcal{N}(t|\boldsymbol{\phi}(x)^\mathrm{T} \boldsymbol{w} , \beta^{-1})
}

とすると,ここからp(t|x,\mathcal{D})

\displaystyle{
p(t|x,\mathcal{D}) = p(t|\boldsymbol{m}^\mathrm{T}\boldsymbol{\phi}(x),\sigma(x)^2) \\
\sigma(x) = \frac{1}{\beta} \bigl( 1 + \boldsymbol{\phi}(x)^\mathrm{T} \boldsymbol{S} \boldsymbol{\phi}(x) \bigr)
}

です.

    1. ビショップ,2006, パターン認識機械学習, 丸善出版, p.p.30 - 31.

モデルエビデンス

 良いモデルとは端的に言ってモデルエビデンス


ME = p(\mathcal{D}|\mathcal{M}_i)

を最大化するような\mathcal{M}_iを選ぶことです.上式は


p(\mathcal{D}|\mathcal{M}_i) = \int p(\mathcal{D}|\boldsymbol{w}, \mathcal{M}_i)p(\boldsymbol{w}|\mathcal{M}_i) d\boldsymbol{w}

と書け,すなわち,事前分布p(\boldsymbol{w}|\mathcal{M}_i)からランダムにサンプリングしたときに\mathcal{D}が生成できる確率を表しています.最大事後確率推定(MAP推定)で推定された係数\boldsymbol{w}_\mathrm{MAP}の周りで確率が鋭く尖っていると仮定すると上式の対数は正規分布では

\displaystyle{
\ln p(\mathcal{D}|\mathcal{M}_i) \approx \ln p(\mathcal{D}|\boldsymbol{w}_\mathrm{MAP}) +\frac{1}{2} \ln \frac{\mathrm{det}\boldsymbol{S}_\mathrm{posterior}}{\mathrm{det}\boldsymbol{S}_\mathrm{prior} }}

と近似できます.左辺第1項は予測分布と観測されたデータの誤差,第2項は事前分布と事後分布の分散の広がり(体積)の比です.実際に変化を見てみましょう

モデルの変化(基底の数)

基底として,xの累乗x_mを3個,4個,10個用意した場合を考えます.テストデータは[0,1]の間にN=20個で,分散\sqrt{1/\beta}=0.2を持っているとします.以下にテストデータと真の分布のグラフを載せます. 20190406231704

ひとまず,事前分布として


p(\boldsymbol{w}) = \mathcal{N}(\boldsymbol{w}|\boldsymbol{0},1/\alpha \beta \boldsymbol{I})

を用いた場合を考えましょう.\alpha = 10^{-8}とします(本来は\betaは未知で,\alpha含めてエビデンス近似などで最適化する必要がありますがここでは割愛します).これによる予測分布を下図に示します.

二次関数までしか使っていないM=3ならともかくM=4,10はそこまで予測として違いがないように思えます.では,\boldsymbol{w}に関する事前分布と事後分布はどうでしょうか?事前分布と事後分布の平均と分散,つまり\boldsymbol{w}_\mathrm{MAP}と周辺化された分散\sigma_m = \sqrt{ \{ \boldsymbol{S} \}_{mm} / \beta}を見てみましょう.

基底の数が多いほど,推定された平均が上下に大きく振動し,事前分布に対して事後分布の分散も対して変わっていないように見えます.もはや近似式が使えるかも怪しいです.モデルエビデンスを計算しましょう.

\displaystyle{
ME = \ln p(\mathcal{D}|\boldsymbol{w}_\mathrm{MAP}) + \frac{1}{2} \ln \frac{\mathrm{det} \beta^{-1} \boldsymbol{S}}{\mathrm{det} (\beta \alpha)^{-1} \boldsymbol{S}_0}  \\
= - \frac{1}{2} \left( N \ln \frac{2 \pi}{\beta} + \beta || \check{\boldsymbol{t}} - \boldsymbol{\Phi} \boldsymbol{w}_\mathrm{MAP} ||^2 \right) + \frac{1}{2} \ln \frac{\mathrm{det} \boldsymbol{S}}{\mathrm{det} \alpha^{-1} \boldsymbol{S}_0}  
}

と計算できます.基底の数ごとのモデルエビデンス(ME)は下のグラフとなります.

20190406231654

M=4を越したあたりから,徐々に減っていっていることがわかります.十分な表現力を持つ関数の数ではM=4が最低なので,当たり前といえば当たり前ですね.一般にモデルエビデンスは基底の数に対して最大値を越すと単調に減少して行く場合が多いです.

モデルの変化(基底の種類)

もし,グラフの分布から,基底の種類を\sin2\pi x1で表されると考えた場合はどうでしょう?回帰の結果と計算されたMEを下に示します. 20190406234753

モデルエビデンスxの累乗の場合よりかなり少なくなっています.このようにテストデータに対して「表現力のいい」関数を選ぶこともモデル選択において大切であることがわかります.

参考文献

    1. ビショップ, 2006, 「パターン認識機械学習 上」,丸善出版,p.p.160 - 164.

追記

4/6:図を修正しました.