「大規模計算時代の統計推論」 を全部Rでやってみる ~第3章~
こんにちは,kur0ckyです.
最近の出社と在宅の比率は半々ぐらいです.緊急事態宣言がでていた最後の方はメンタルきつかったのですが,今は落ち着いています.あまりの無気力さに,スーパーに行っても何も買う気が起きなかったりしていました.笑
第3章は,「ベイズ派的な推論」です. 前章は頻度派のイントロでしたが,今回はそのベイズ版です. 今回も前回に引き続き解析は少ないですが,それ以外にも初学では辛そうな知識や抽象的な記述が多いところは備忘録的に整理しています. 誤りが多分にあると思うので,見つけたときは石を投げてください. 現在読むだけなら半分ぐらいきてるのですが,この本はすごく勉強になってるので,買うと良いと思います.
- 作者:エフロン,ブラッドレイ,ヘイスティ,トレバー
- 発売日: 2020/07/30
- メディア: 単行本
前回記事はこちらです.
それではやっていきます.
本章は次のような記述からはじまります.
ベイズ派の推論法は頻度主義と真っ向から対立している.というのは言いすぎだとしても,少なくとも直交はしている.
(中略)
2つの哲学の美点を合体させるための苦労は,より込み入ったデータ集合を扱わなければならない現代では,より重大になっている.
いきなり読み進めるのが楽しみになるような文言ですね. 「直交する」というのはどのようなことなのか,いろいろ考えてしまいますが,進めていきます.
はじめに
頻度派においてもベイズ派においても,確率密度の族がその根幹にあります. まずはその数式の準備がされています.
標本空間,パラメータ空間,密度関数からなる族をとします. はそれぞれ空間内のひとつの値です. この族に加え,ベイズ派ではさらに事前知識の仮定をおきます. つまりです. そもそも「事前知識」とはなにか,という率直な疑問についてはさほど触れられておらず,核となるベイズの定理の説明に進みます. どのような場合に事前知識が明らかかは,この章のいくつかの事例を通して触れられています.
ベイズの定理は次式で与えられます.
ベイズの定理や事後分布についての説明は,どの基本的な統計の本にも載っているような事項なので,この記事では割愛します. 上式より明らかですが,事後分布では観測値に固定されています.一方で,はに渡って変化する変数として扱っています. これは頻度派とはちょうど逆になっていますね. さらに,は固定されているので,分母のは,定数扱いすることができます. 事後密度の面積を1にするための正規化定数と捉えることもできます. また,任意の異なるパラメータに対しては共通しているため,それぞれのもとでの事後密度の比は次のようになります.
つまり,を考えずに両者を比較・評価することができます. この便利な形(事後オッズ比)は,ベイズ流の検定やMCMCのアルゴリズムなどでしばしば利用されている気がします.
3.1 2つの例
まず3.1節では,ベイズ則がどのように使われるか,2つの例を示してその便利さが示されています.
物理学者の双子の例
まず1つ目は,ある物理学者が双子の男の子を身篭っていると知ったとき,その双子が二卵性双生児であるか一卵性双生児であるかを推論する問題です. 簡潔に整理してみます.
一卵性では,基本的に同性になります(少し調べてみると,稀に違うときもあるようです). 一方で二卵性では,同性になる確率も異性になる確率も1/2になります. 一卵性になるか二卵性になるかのベルヌーイ試行で,パラメータが異なるイメージです. また,双子が同性であるかどうかも二者択一であるため,こちらもベルヌーイ試行にしたがいます. そこで,同性になるとき1,異性になるとき0になるような確率変数を考えます.
また,上述の事実より,双子が同性になるかどうかの確率を支配する母数は,双子が一卵性か二卵性かによって異なります. 一卵性のときを,二卵性のときをとしましょう. すると,この問題は,双子が同性である,つまりを観測したときに,の事後確率はそれぞれどれほどになるか,という問題で考えることができます. つまり,次の二つを比較すれば良いことになります.
- .
- .
これらは,が自然の摂理から明らかであること,であることなどを利用すれば求めることができます. 本書では前述の事後オッズ比を求めており,「双子の男の子」という観測のもとで両者の事後確率が等しいことを示しています.
この例を整理した図3.1は解析ではないためこの記事で再現はしませんが,これはの同時分布を整理したものです. もも2つの値しか取らないため,このような図で整理することができますが,次元が大きくなったり連続変数になったりするとこのような整理は不可能です. しかし,この図が分かり易くしているような,同時確率と周辺確率,事後確率,事前確率をいったりきたりすることはベイズ派,ひいては頻度派でも超重要と思われます. 私が技術書や論文を読むときは数式をさくっと追うだけになってしまいがちなのですが,逃げることなくしっかりやりたいです(自戒).
テストの点数の例
2つ目の例は,22名の学生に対して実施した2科目のテスト点数の相関係数についての例です. 大きさ22の標本から普通にもとめた相関係数は0.498になりますが,ありとあらゆる学生に対して実施したときの相関係数を知りたい,というのは自然な欲求です.
これは後の5章で議論されることですが,2つの科目の点数の同時分布が2変量正規分布に従うとしたとき,その相関係数の分布は次式になります.にはデータから推定した相関係数の値を差し込み (plug-in) ます.
この例では,事前知識の入れ方として,次の3つの場合を紹介しています.
- 無情報事前分布
- .
- 三角事前分布
- .
- ジェフリーズの事前分布
- .
式より,無情報事前分布はがにわたり一様に分布している,平坦な分布です. 三角事前分布は,0を中心として三角形の形をしている分布です. ジェフリーズの事前分布に関しては知識がないと何が何やら分からないと思うので.次で整理したいと思います.
それでは,事後分布はで計算できることを思い出し,それぞれの事前分布から計算してみましょう. ここで,は正規化定数とみなせるため,あとで積分して1に調整してしまうことにしていったん放置しましょう. それぞれ計算して,各事前分布を置いたときの事後確率は次のように表せます. には22を代入してあります.
無情報事前分布の場合
三角事前分布の場合
ジェフリーズの事前分布の場合
図3.2
以上,事後分布がもとまったので,Rで作図していきましょう.
tidyverse
を読み込み,データをダウンロードします.
また,差し込む相関係数の推定値も算出しておきます
library(tidyverse) scores <- read_delim("https://web.stanford.edu/~hastie/CASI_files/DATA/student_score.txt", delim = " ") # 相関係数の算出 theta_hat <- cor(scores$mech, scores$vecs)
の積分の部分はすべてに共通しているため,f_denom()
としてわけて考えます.
残りの部分はそれぞれf_***_nume()
として定義します.
ネーミングは雑に「分子」「分母」から取りましたw.
f_denom()
では数値積分を行う必要があるのですが,ここではひと工夫必要です.
Rで数値積分するにはintegrate()
を使えば良いのですが,を与えつつの関数として扱いたいため,さらに関数を戻り値にしています.
# 共通するwで積分する部分 f_denom <- function(theta){ return(function(w){1 / (pi * (cosh(w) - theta*theta_hat)^21) }) } # 無情報事前分布 f_flat_nume <- function(theta){ 10*(1-theta^2)^10.5*(1-theta_hat^2)^9 } f_flat <- function(theta){ f_flat_nume(theta) * integrate(f_denom(theta), 0, Inf)$value } # 三角事前分布 f_triangle_nume <- function(theta){ 20*(1-abs(theta))*(1-theta^2)^10.5*(1-theta_hat^2)^9 } f_triangle <- function(theta){ f_triangle_nume(theta) * integrate(f_denom(theta), 0, Inf)$value } # ジェフリーズの事前分布 f_jeff_nume <- function(theta){ 20*(1-theta^2)^9.5*(1-theta_hat^2)^9 } f_jeff <- function(theta){ f_jeff_nume(theta) * integrate(f_denom(theta), 0, Inf)$value }
次に,放置していた正規化定数を計算しましょう.
確率密度は面積が1にならなくてはならないので,最後に正規化定数で割らなければなりません.
関数を入れ子にしているためか,integrate()
ではうまくいかなかったので,重積分も可能なcubature::adaptIntegrate()
を使います.
必要に応じてインストールinstall.packages("cubature")
します.
library(cubature) norm_const_flat <- adaptIntegrate(f_flat, -1, 1)$integral norm_const_triangle <- adaptIntegrate(f_triangle, -1, 1)$integral norm_const_jeff <- adaptIntegrate(f_jeff, -1, 1)$integral
それでは作図します.ggplot2::stat_function()
を使おうかとも思ったのですが,面倒だったので作った関数で値を計算して素直にgeom_line()
で書きました.
また,軸ラベルに数式を書くところではTeX記法を使いたかったので,latex2exp::TeX()
を使っています.expression()
で頑張っても良いでしょう.
これも必要に応じてインストールしてください.
library(latex2exp) tibble(theta = seq(-1, 1, 0.01)) %>% mutate(p_flat = map_dbl(theta, f_flat) / norm_const_flat, p_triangle = map_dbl(theta, f_triangle) / norm_const_triangle, p_jeff = map_dbl(theta, f_jeff) / norm_const_jeff) %>% gather(key, p, -theta) %>% ggplot(aes(theta, p, group = key))+ geom_line(aes(linetype = key))+ geom_vline(xintercept = theta_hat)+ scale_x_continuous(limits = c(-0.3, 1))+ labs(x = TeX('$\\theta$'), y = TeX('$g(\\theta | \\hat{\\theta})$'), linetype = "事前分布")
無事綺麗な図を書くことができました.
3.2 無情報事前分布
さて,テストの例では,特に2科目の点数の相関について事前知識があるわけではありませんでした. このように,役立つ事前情報に恵まれない場合でも,ベイズ即を使いたいときがあるかもしれません. その一つの方法が無情報事前分布(uninformative prior)です.
無情報事前分布は前節でも使いましたが,事前知識を持たないことを再現したいので,いたるところで一様な分布ということになります. ちなみに,これは多くの人が聞いたことがあると思いますが,無情報事前分布を用いた場合のMAP推定値は,次式のように最尤推定値に等しくなります. 読み進めるにつれまた扱う機会があると思います.
しかし,無情報事前分布には一つ欠陥があります. パラメータ空間がであるとき,無情報事前分布を一様分布としてしまうとになってしまうのです. これは確率の公理を満たさず,improper prior distributionとよばれます.
十分広い閉区間に局所一様分布を考えれば問題ないですが,これも変数変換などによって一様性が失われてしまいます. これを解決するための手法に,先ほどもでてきた「ジェフリーズの事前分布」があります.
ジェフリーズの事前分布
初見では理解できないと思ったので簡単に解説を入れることにします.
ジェフリーズの事前分布とは,特定の事前分布のことを指すのではなく,事前分布を構成するための「規則」です. そして,この「規則」とは,「パラメータを変数変換しても事後分布への影響が不変になるように事前分布を選択すること」です.
どうしてそのようになるかの証明は後回しにし,まずはジェフリーズの事前分布の定義を紹介します.
ここで,はのフィッシャー情報量です. つまり,ジェフリーズによる「規則」は,「パラメータのフィッシャー情報量に比例するように事前分布を選ぶ」という規則でもあります. フィッシャー情報量は4章ででてくるのでここでは割愛しますが,簡単に言うと「対数尤度の偏微分(スコア関数)の分散」です. 推定量の良さを評価する一手法クラメル・ラオの下限に使われることで知っている人も多いと思います.
詳しく知りたい方は以下の記事などを読むとわかりやすい気がします.
定義を紹介したところで,本当にパラメータの変数変換に対して不変であるのか,示したいと思います.
上の式変形より,から変数変換をしてもがいえたので,変数変換に対して不変であるとわかりました. 半ばトップダウン的な形で示しましたが,なんだか美しいですね. フィッシャー情報量の平方根に比例さえするように事前分布を選べば,変数変換しても問題ないというわけです. パラメータが複数の場合は行列を考えなければいけないのですが,ここでは割愛します.
3.3 頻度派的な推論の欠陥
この節では,3つの例を通して,頻度派的推論の弱点を指摘しています.
検診員の例
この例に関しては正直よくわかりませんでしたが,なんとか次のように考えました. 間違っていたら教えて欲しいです.
まずは次の記述についてです.
次の日,彼は電圧計に不具合を見つけた.電圧が100を超えた場合はすべてと報告されたかもしれないというのである.頻度派統計学者風に考えると,は真の期待値に対してもはや不変ではない.
は,計測値がになるように調整されている,という前提のもとで不偏な推定量であるため,その前提が覆された以上不偏な推定は行えない,ということでしょうか. ここまでは何となく良さそうな気がします.
つづいて次の記述があります.
ベイズ派の統計学者風に考えてみる. 任意の事前密度に対して,12個の計測値からなるベクトルにおける事後密度は,実際に観測されたデータのみに依存し,「観測されたかもしれない他の潜在的なデータ集合」には依存しない.
これについてはよく分かりませんでした. 計測値が正規分布に従うという前提がある以上,尤度の計算には正規分布の密度関数を使うことになると思います. その前提が崩れた以上,事後密度が「データのみに依存する」からといって,を事後期待値として報告することには依然として問題がある気がしました.
途中で実験をやめる例
典型的なp-hackingの例です. この例では,実験者が予定していた30回の実験を待たずに,「20回目であわよく有意な結果が得られていれば実験をそこで終わりにしよう」と結果を盗み見てしまいました. そのため,結局30回目で有意になり棄却できたとしても,もはやこの実験の有意水準は元の値ではなくなってしまった,ということになります. はじめから,「20回もしくは30回目のいずれかで有意になったときに実験を停止する」という規則で検定を設計しなければなりません. p-hackingはよくないですね. 結果が有意になるまでデータを採り続ける,というのも同様です.
頻度派的な推論ではついこのようなことが出来てしまいますが,ベイズ派では実験が早く中止になったかどうかに関わらず,事後分布を見ればよいため,寛容である.との主張がなされています..
この例はランダムネスのみに依存して再現不可能なので,割愛します.
prostateデータの例
この例では,前立腺がんの研究に関するデータを扱います. 52人の患者と50人の健常者に対し,6033個の遺伝子について値を測定したデータです. 解析例では,すべての遺伝子について,患者と健常者で平均の差の検定を行っています. また,その検定統計量をなぜか標準偏差1の正規分布に従うように変換しています. 理由はよくわかっていないのですが,Efron(2010, Section 7.4) を読めば分かるようです. ここでは深く追いません.
https://www.cambridge.org/core/books/largescale-inference/A0B183B0080A92966497F12CE5D12589
まずはデータを読み込みましょう.
prostmat <- read_csv("https://web.stanford.edu/~hastie/CASI_files/DATA/prostmat.csv") prostz <- scan("https://web.stanford.edu/~hastie/CASI_files/DATA/prostz.txt", what = double())
それでは解析していきます.
すべての遺伝子に対して検定をしなければならないので,遺伝子ごとにインデックスgene_index
をつけ,nest()
し,map()
で回すのが早そうですね.
また,検定統計量を正規分布に従うように直す作業はqnorm(pt(statistic, df = 100)))
でやっています.
pt()
で自由度100の分布における下側確率を求めた後,qnorm()
で,その下側確率に対応する正規分布での値を求めています.
prostate_test <- prostmat %>% mutate(gene_index = 1:n()) %>% gather(key, value, -gene_index) %>% mutate(key = str_remove_all(key, "_.*")) %>% group_by(gene_index, key) %>% nest() %>% ungroup() %>% mutate(data = map(data, ~ unname(unlist(.x)))) %>% spread(key, data) %>% mutate(t_test = map2(cancer, control, ~ t.test(.x, .y, paired = F, var.equal = F)), statistic = map_dbl(t_test, ~ .x$statistic), x = qnorm(pt(statistic, df = 100))) print(prostate_test)
出力すると次のようになります.
私はnest()
とmap()
信者なのですが,このように管理できるのはたとえば交差検証やtarget encodingでもめちゃくちゃ便利です.
> print(prostate_test) # A tibble: 6,033 x 6 gene_index cancer control t_test statistic x <int> <list> <list> <list> <dbl> <dbl> 1 1 <dbl [52]> <dbl [50]> <htest> 1.48 1.47 2 2 <dbl [52]> <dbl [50]> <htest> 3.70 3.57 3 3 <dbl [52]> <dbl [50]> <htest> -0.0278 -0.0278 4 4 <dbl [52]> <dbl [50]> <htest> -1.14 -1.13 5 5 <dbl [52]> <dbl [50]> <htest> -0.141 -0.140 6 6 <dbl [52]> <dbl [50]> <htest> 0.963 0.959 7 7 <dbl [52]> <dbl [50]> <htest> 1.07 1.07 8 8 <dbl [52]> <dbl [50]> <htest> -1.30 -1.29 9 9 <dbl [52]> <dbl [50]> <htest> -1.24 -1.23 10 10 <dbl [52]> <dbl [50]> <htest> 1.18 1.18 # … with 6,023 more rows
それでは作図してしまいます.
prostate_test %>% ggplot(aes(x))+ geom_histogram(colour = "gray25", fill = "green", bins = 50)+ labs(x = "効果量の推定値", y = "度数")
本書では,遺伝子610が示す最も大きい値5.29をとりあげ,これは過大評価されており,下方修正すべきとしています. 「最大値を取る遺伝子をとってくる」という作業をしているので,これは順序統計量と比較しなければならないということでしょうか. この例は15章でまた出てくるらしいので,そのときを楽しみに待つことにします. ちなみに私はまだそこまで読み進められていません笑
3.4 ベイズ派と頻度派の対照表
この節では,頻度派とベイズ派の比較が行われています. 3章まででは具体的な手法がそれほど解説されているわけではないので,抽象的な表現もあります.むしろ読み進めた後に立ち返ると良いのでは,と個人的には思いました. 私の実力不足もあり完全には噛み砕けていませんが,ここでは備忘録的に少し整理します. 解析はありません.
両主義の比較は,次の文章を起点とし,細かいスタンスの違いや良し悪しを列挙する形で行われていきます.
どちらも式(3.1)の確率分布族から出発しており,いわば最初は同じ場所で仲良く遊んでいたのだが,そのうちに直交した方向に歩み始める.
(中略)
ベイズ推論は垂直方向に,すなわちを固定した上で事後分布に従って進むが,頻度派の考え方は水平方向である.すなわち,を固定した上で,を変化させて考える.
ここまで整理してきても何が「直交」なのかいまいち分かってないのですが,結局単純に次のような整理を比較させただけでしょうか.
- ベイズ派:得られたデータを所与とし,確率モデルと事前分布の仮定のみで事後分布を推定する.
- 頻度派:データの背後に一定の確率分布が存在すると仮定し,その着目する性質と,それを推定する想定するアルゴリズムで議論する.議論が十分になされた後,選択したアルゴリズムにデータを代入する.
そこからの比較は箇条書き的に行われており,超要約すると次のような感じです.
事前分布が具体的に想定される場合は躊躇せずベイズの定理を使える.そうでない場合はジェフリーズの事前分布のようなテクニックもあるが,論理的一貫性の魅力(???)は享受できない.
良さげな事前分布がない場合,どうしても恣意性がつきまとう.頻度派は常に客観性を要求する.
ベイズ派では,一度事前分布を定めてしまえば,(MCMCや変分推論のようなアルゴリズムはさておき)事後分布だけに着目すればよい.頻度派で選択するアルゴリズムは全てそこから計算できる.
一方で頻度派では,頻度派では,事前分布を選ぶ代わりに,注目したい母集団の性質ごとにアルゴリズムを選択し,その確率的性質に関して議論する.異なる問題に対しては異なるアルゴリズムを考える.想定される分布族の元でアルゴリズムの性質を議論し,比較検討しなければならない.そのため,分布とアルゴリズムに対しone-to-oneに設計されるため,場当たり的になる.手始めに無情報事前分布をおいて行うベイズ推論に人気があるのも,ここを避けれることが大きな理由.
事前分布の選択こそベイズの課題.頻度主義はより慎重な姿勢.母数がどんなものであれ破綻しないように議論する.もしくは特定の状況下であればうまくいく,というような議論が行われる.
動的な更新が必要な状況下ではベイズは魅力的.物事が進展するにつれ確率が更新されていくのは自然だし,共役事前分布などが大活躍する.
しばしば出てくるベイズの「論理的一貫性」という表現について,あまりピンときていません. 単に,事前分布と確率モデルから事後分布の推論のみの議論だけで最後までいける,ということでしょうか.
おわりに
3章も無事まとめ終わりましたが,割と広い知識が求められているので良い復習になっています. 4章以降は具体的な論理が展開されていくので,はやく次の記事を出せるようにしたいです. 割と時間がかかってしまっています... 社会人になり,どのように時間を使っていくか,悩ましい限りです.