ベイズ統計_03_尤度

ベイズ統計 03 尤度 分寺 杏介 神戸大学大学院 経営学研究科  bunji@bear.kobe-u.ac.jp ※本スライドは，クリエイティブ・コモンズ 表示-非営利 4.0 国際 ライセンス（CC BY-NC 4.0）に従って利用が可能です。

• https://creativecommons.org/licenses/by-nc/4.0/

前回のおさらい ベイズの定理はさすがにもう覚えましたか？ A) ベイズの定理の式を書いてみてください。 B) ベイズの定理に登場する４つの確率分布について説明してください。 C) 確率分布の「正規化定数」と「カーネル」を説明してください。 ※次ページ以降に解答例が載っています。自力で解答したい人は進まないように！！ 03 尤度 2

ベイズの定理と確率分布 ここからはデータを𝑌 ，パラメータを𝜃 とする ベイズの定理 𝑃 𝑌𝜃 𝑃 𝜃 𝑃 𝜃𝑌 = 𝑃 𝑌 【事後分布】 Posterior Distribution 最終的に求めたいやつ。 あるデータが所与のときに，求めたいパラメータの確率を表す。 03 尤度 3

ベイズの定理と確率分布 ここからはデータを𝑌 ，パラメータを𝜃 とする ベイズの定理 𝑃 𝑌𝜃 𝑃 𝜃 𝑃 𝜃𝑌 = 𝑃 𝑌 【事前分布】 Prior Distribution データが与えられる前から持っている，パラメータ自体の確率分布。 なおベイズ統計では主観確率分布を置くことが可能。つまり • 去年の国勢調査から得られた分布，みたいなデータに基づくものでも • 「俺はこう思う」というなんの根拠のないものも（理論上）許してくれる。 03 尤度 4

ベイズの定理と確率分布 ここからはデータを𝑌 ，パラメータを𝜃 とする ベイズの定理 𝑃 𝑌𝜃 𝑃 𝜃 𝑃 𝜃𝑌 = 𝑃 𝑌 【尤度】 Likelihood 今回はこいつのお話（おさらい）をしていきます。 03 尤度 5

ベイズの定理と確率分布 ここからはデータを𝑌 ，パラメータを𝜃 とする ベイズの定理 𝑃 𝑌𝜃 𝑃 𝜃 𝑃 𝜃𝑌 = 𝑃 𝑌 【周辺尤度・エビデンス】 Marginal Likelihood こうやって見ると「データが得られる確率」ということになる。 とりあえず今はあまり気にしなくていい。来る時が来たらその意味がわかるだろう… 03 尤度 6

何が正規化定数で何がカーネルか 𝑃 𝑥|𝜇, 𝜎 = = 正規化定数 × 1 × 2𝜋𝜎 2 𝑥の値によらず一定な部分 ベイズの定理 カーネル 𝑥−𝜇 2 exp − 2𝜎 2 𝑥の値によってかわる部分 𝑃 𝜃𝑌 ∝𝑃 𝑌𝜃 𝑃 𝜃 𝑃 𝑌𝜃 𝑃 𝜃 𝑃 𝜃𝑌 = 𝑃 𝑌 実際に推定する際には カーネルのみを考えることが多い 03 尤度 7

１ 尤度とは 03 尤度 8

条件付き確率 𝑃 𝑌 𝜃 は「パラメータが𝜃のとき，データ𝑌が得られる確率」 【尤度】 例 当たり確率 𝜋 = 0.4のクジを3回ひいたとき，ちょうど 𝑘 回当たる確率は？ ▶ 二項分布に従うとすると… 𝑃 𝑘 𝜋, 𝑛 = 𝑛𝐶𝑘 𝜋 𝑘 1 − 𝜋 𝑛−𝑘 𝑌, 𝜃 = (𝑘, 𝜋) 今回はこいつのお話（おさらい）を ※ 今回 𝑛 は所与なのでパラメータではない 確率の とベイズの定理 𝜋 = 0.4に固定された状態で 𝑘の値を動かすと確率がどう変化するか 【ちょっとだけ練習】 𝜋 = 0.4のときに3回中1回当たる確率は？ A 03 尤度 およそ25.1% 9

問題設定を変えてみると 例 クジを3回ひいたらちょうど1回当たった。このときクジの当たり確率 𝜋 は？ ▌問題を簡略化して考えてみましょう • クジの当たり確率には3種類の 設定 があり毎日異なっている。 • 以下のうち，今日の 設定 はどれでしょうか？ 設定1 当たり確率は0.2 ▶ 二項分布 𝑃 𝑘 𝜋 = 0.2, 𝑛 = 3 設定2 当たり確率は0.5 ▶ 二項分布 𝑃 𝑘 𝜋 = 0.5, 𝑛 = 3 03 尤度 設定3 当たり確率は0.8 ▶ 二項分布 𝑃 𝑘 𝜋 = 0.8, 𝑛 = 3 10

各設定のもとで「データが得られる確率」を計算 • 以下のうち，今日の 設定 はどれでしょうか？ 設定1 設定2 設定3 当たり確率は0.2 ▶ 二項分布 𝑃 𝑘 𝜋 = 0.2, 𝑛 = 10 当たり確率は0.5 ▶ 二項分布 𝑃 𝑘 𝜋 = 0.5, 𝑛 = 10 当たり確率は0.8 ▶ 二項分布 𝑃 𝑘 𝜋 = 0.8, 𝑛 = 10 1 2 3C1 0.2 0.8 = 0.384 1 2 3C1 0.5 0.5 = 0.375 1 2 3C1 0.8 0.2 = 0.096 この3つの中であれば 設定1 が最も 𝑘 = 1 を生成しうる設定である！ 03 尤度 11

連続変数に拡張してみると 例 クジを3回ひいたらちょうど1回当たった。このときクジの当たり確率 𝜋 は？ ▌問題をより一般化してみると • クジの当たり確率は0から1のうちのいずれかの値 連続変数であればいくらでも細かい値を考えることができる • いわば，無限通りの 設定 があるような状態 𝜋 = 0.001, 𝜋 = 0.27482763, 𝜋 = 0.5012837 … ▌ 先ほどと同様に二項分布に従って𝑃 𝑘 = 1 𝜋, 𝑛 = 3 をすべての 𝜋 で計算すると 𝑃 𝑘 = 1 𝜋 = 0.2 = 0.384 𝑃 𝑘 = 1 𝜋 = 0.5 = 0.375 データ (𝑘 = 1) を所与としたときの パラメータ 𝜋 の関数とみなせる 𝑃 𝑘 = 1 𝜋 = 0.8 = 0.096 条件付ける側と付けられる側を ひっくり返した関数 𝑃 𝜋 𝑘, 𝑛 について考えていきます 03 尤度 12

二項分布の場合 𝑃 𝑘𝜋 𝑃 𝜋 = もし𝜋の値を一つに決めたら， 𝑃 𝑘, 𝜋 = 𝑃 𝑘 𝜋 𝑛 は所与のため省略しています 𝑃 𝑘, 𝜋 ※もし𝑘, 𝜋が両方とも 固定されれば，左右の 2つは全く同じになる = 𝑃 𝜋 𝑘 𝑃(𝑘) もし𝑘が所与だったら， 𝑃 𝑘, 𝜋 = 𝑃 𝜋 𝑘 𝜋がある値のときの𝑘の確率 𝑘がある値のときの𝜋の確率 𝜋をある値に固定した状態で 𝑘をある値に固定した状態で 𝑘を動かすと確率がどう変化するか 𝜋を動かすと確率がどう変化するか 標本理論では母数はある固定の値 ベイズでは母数は変動するもの(分布) データはサンプリングで変動するもの データは所与のもの 見方が違うだけで，同じ確率𝑃 𝑘, 𝜋 に基づいた話をしている 03 尤度 13

一般化すると 𝑃 𝑌𝜃 𝑃 𝜃 = もし𝜃の値を一つに決めたら， 𝑃 𝑌, 𝜃 = 𝑃 𝑌 𝜃 𝑃 𝑌, 𝜃 ※もし𝑌, 𝜃が両方とも 固定されれば，左右の 2つは全く同じになる = 𝑃 𝜃 𝑌 𝑃(𝑌) もし𝑌が所与だったら， 𝑃 𝑌, 𝜃 = 𝑃 𝜃 𝑌 𝜃がある値のときの𝑌の確率 𝑌がある値のときの𝜃の確率 𝜃をある値に固定した状態で 𝑌をある値に固定した状態で 𝑌を動かすと確率がどう変化するか 𝜃を動かすと確率がどう変化するか 標本理論では母数はある固定の値 ベイズでは母数は変動するもの(分布) データはサンプリングで変動するもの データは所与のもの 見方が違うだけで，同じ確率𝑃 𝑌, 𝜃 に基づいた話をしている 03 尤度 14

確率（密度）関数を別の見方で パラメータが所与の場合（確率関数） ▌𝑛 = 3 の二項分布の場合 𝑃 𝑌 = 𝑘 𝜋 = 3𝐶𝑘 𝜋 𝑘 1 − 𝜋 3−𝑘 = 𝑃 𝜃 = 𝜋 𝑘 データが所与の場合（尤度関数） 𝑃(𝑘 = 0) 𝑃(𝑘 = 1) 𝑃(𝑘 = 2) 𝑃(𝑘 = 3) 𝜋 = 0.0 1.000 0.000 0.000 0.000 𝜋 = 0.1 0.729 0.243 0.027 0.001 𝜋 = 0.2 0.512 0.384 0.096 0.008 𝜋 = 0.3 0.343 0.441 0.189 0.027 𝜋 = 0.4 0.216 0.432 0.288 0.064 𝜋 = 0.5 0.125 0.375 0.375 0.125 𝜋 = 0.6 0.064 0.288 0.432 0.216 𝜋 = 0.7 0.027 0.189 0.441 0.343 𝜋 = 0.8 0.008 0.096 0.384 0.512 𝜋 = 0.9 0.001 0.027 0.243 0.729 𝜋 = 1.0 0.000 0.000 0.000 1.000 03 尤度 確率関数 𝑃 𝑌=𝑘𝜋 尤度関数 𝑃 𝜃=𝜋𝑘 15

ただし「確率」ではない パラメータが所与の場合（確率関数） ▌𝑛 = 3 の二項分布の場合 𝑃 𝑌 = 𝑘 𝜋 = 3𝐶𝑘 𝜋 𝑘 1 − 𝜋 3−𝑘 = 𝑃 𝜃 = 𝜋 𝑘 データが所与の場合（尤度関数） 𝑃(𝑘 = 0) 𝑃(𝑘 = 1) 𝑃(𝑘 = 2) 𝑃(𝑘 = 3) （確率関数なので言うまでもなく） 1.000 0.000 0.000 0.000 𝜋 = 0.0 𝜋 が所与のもとでは 0.729 0.243 0.027 0.001 𝜋 = 0.1 𝜋 = 0.2 0.512 0.384 0.096 0.008 𝑛 𝜋 = 0.3 0.343 0.441 0.189 0.027 ෍𝑃 𝑌 = 𝑘 𝜋 = 1 𝜋 = 0.4 0.216 0.432 0.288 0.064 𝑘=1 𝜋 = 0.5 0.125 0.375 0.375 0.125 𝜋 = 0.6 0.064 0.288 0.432 0.216 𝜋 = 0.7 0.027 0.189 0.441 0.343 න 𝑃 𝜃 = 𝜋 𝑘 = 1 とは限らない 𝜋 = 0.8 0.008 0.096 0.384 0.512 𝜋=0 𝜋 = 0.9 0.001 0.027 0.243 0.729 𝜋 = 1.0 0.000 0.000 0.000 1.000 03 尤度 𝑘 が所与のもとでは 1 𝑃 𝜃 = 𝜋 𝑘 は確率とは違う！ 𝐿 これを尤度 (likelihood) と呼びます 16

尤度関数の使い方｜頻度主義の場合 ▌尤度関数は「所与のデータ𝑌は，どの 設定 から発生しやすいか」を表すもの 頻度主義的には，パラメータの「真の値」＝「真の 設定 」がただ一つ存在する ▶ 統計的推測としては，最も発生しやすい 設定 を選ぶのが良かろう ▌最尤推定 (maximum likelihood estimation: MLE) 例 クジを3回ひいたらちょうど1回当たった。このときクジの当たり確率 𝜋 の最尤推定値は？ 尤度関数は𝐿 𝜋 𝑘 = 1 = 3𝐶1 𝜋 1 1 − 𝜋 2 = 3𝜋 1 − 𝜋 2 【左図】 1 3 この尤度関数は， 𝜋 = のときに最大値を取るので 最尤推定値は 𝜋ො = 1 3 「３回中１回当たった」という結果は 𝜋ො = 1 の設定から発生したと考えるのが 3 一番尤もらしい推論ですねぇ 03 尤度 17

尤度関数の使い方｜ベイズ統計の場合 ▌尤度関数は「所与のデータ𝑌は，どの 設定 から発生しやすいか」を表すもの ベイズ統計的には，各 設定 に対する信念の強さを確率分布として表す ▶ ベイズ統計では，尤度関数がそのまま各 設定 に対する信念の更新式となる 例 クジを3回ひいたらちょうど1回当たった。このときクジの当たり確率 𝜋 の事後確率は？ 尤度関数は𝐿ベイズの定理と確率分布 𝜋 𝑘 = 1 = 3𝐶1 𝜋 1 1 − 𝜋 2 = 3𝜋 1 − 𝜋 2 【左図】 ここからはデータを ，パラメータを とする これがそのまま 「データによる信念の更新」 を表す 03 尤度 【尤度】 18

前回登場した面積による考え方 ▌「すべてのあり得る 設定 から，どれが最もありえそうか」を判断 設定1 設定2 設定3 𝜋 = 0.2 𝜋 = 0.5 𝜋 = 0.8 【𝑘 = 1 になるパターンは】 得各 ら設 れ定 るの 確も 率と （で 尤そ 度の ）デ ー タ が 𝑘=0 𝑘=0 𝑘=1 (0.375) 𝑘=1 (0.384) 𝑘=2 𝑘=2 𝑘=3 𝑘 = 1 (0.096) 設定1 𝜋 = 0.2のとき 𝑘=1 (0.384) 設定2 𝜋 = 0.5のとき 𝑘=1 (0.375) 設定3 𝜋 = 0.8のとき 𝑘 = 1 (0.096) 𝑘=2 𝑘=3 03 尤度 19

前回登場した面積による考え方 ▌「すべてのあり得る 設定 から，どれが最もありえそうか」を判断 設定1 設定2 設定3 𝜋 = 0.2 𝜋 = 0.5 𝜋 = 0.8 【実際に 設定1 なのは】 得各 ら設 れ定 るの 確も 率と （で 尤そ 度の ）デ ー タ が 𝑘=0 𝑘=0 𝑘=1 (0.375) 𝑘 = 1 (0.096) 𝑘=1 (0.384) 𝑘=2 ＝0.449 𝑘=1 (0.384) 𝑘=2 𝑘=2 𝑘=3 𝑘=1 (0.384) 𝑘=3 03 尤度 ＋ 𝑘=1 (0.375) ＋ 𝑘 = 1 (0.096) 20

事前の信念による調整もあり ▌（例）過去の傾向から，今日は 設定1 𝑃 𝜋 = 0.2 = 0.25 設定2 𝑃 𝜋 = 0.5 = 0.5 設定2 である可能性が高いとしたら… 設定3 𝑃 𝜋 = 0.8 = 0.25 【実際に 設定1 なのは】 得各 ら設 れ定 るの 確も 率と （で 尤そ 度の ）デ ー タ が 𝑘=0 𝑘=0 𝑘=1 0.25×0.096 𝑘=1 𝑘=1 0.5×0.375 0.25 × 0.384 𝑘=2 𝑘=1 𝑘=1 0.25 × 0.384 𝑘=2 𝑘=2 𝑘=3 0.25 × 0.384 𝑘=3 03 尤度 ＋ 𝑘=1 0.5×0.375 ＝0.312 ＋ 𝑘=1 0.25×0.096 21

事後分布は「確率分布」なので… ベイズの定理 𝑃 𝑌𝜃 𝑃 𝜃 𝑃 𝜃𝑌 = 𝑃 𝑌 【実際に 設定1 なのは】 𝑘=1 0.25 × 0.384 そのうち特定の可能性が占める割合 𝑃 𝑌𝜃 𝑃 𝜃 𝑃 𝑌 = න 𝑃 𝑌 𝜃 𝑃(𝜃) この操作によって 事後分布は 確率分布になる ＝ 𝑘=1 0.25 × 0.384 Θ 𝑘 = 1 というデータが得られる すべての可能性の和 03 尤度 ＋ 𝑘=1 0.5×0.375 ＝0.312 ＋ 𝑘=1 0.25×0.096 22

ベイズ統計と最尤法 ということでベイズ推定は𝑃 𝜃 𝑌 ∝ 𝑃 𝑌 𝜃 𝑃 𝜃 が基本 正規化定数 𝑃(𝑌) を無視しても 事後分布の形だけは求められる ▌ですが，もし事前分布に一様分布を設定すると… 一様分布は𝜃の値によらず一定な分布 この部分がなくなっても分布の形は同じになる 𝑃 𝜃 𝑌 ∝ 𝑃 𝑌 𝜃 になる 𝑃 𝜃 𝑌 ∝ 𝐿 𝜃 𝑌 と書き換えられる この場合，尤度関数がそのまま事後分布になる 03 尤度 もしも事後分布の最大値を 点推定値として利用する場合には 最尤推定と同じ結果になる，ということ 23

複数のデータがあるときの尤度関数 例 5人がクジを3回ずつひいたところ，それぞれ(1,2,2,1,0)回当たった。 このときクジの当たり確率 𝜋 の事後確率は？ もしサンプリングが独立であれば 例えばクジの総数が決まっている場合，前の人の結果によって当たり確率が変動する 独立ではないケース 5 𝐿 𝜋 𝒀, 𝑛 = 𝑃 𝑌1 𝜋, 𝑛 × ⋯ × 𝑃 𝑌5 𝜋, 𝑛 = ෑ 𝑃 𝑌𝑖 𝜋, 𝑛 𝑃 𝑌𝑖 𝜋, 𝑛 = 𝑛𝐶𝑌𝑖 𝑖=1 𝜋 𝑌𝑖 1 − 𝜋 𝑛−𝑌𝑖 【全部書き下すと】 𝐿 𝜋 𝒀, 𝑛 = 3𝐶1 𝜋1 1 − 𝜋 2 × 3𝐶2 𝜋 2 1 − 𝜋 1 × 3𝐶2 𝜋 2 1 − 𝜋 1 × 3𝐶1 𝜋 1 1 − 𝜋 2 × 3𝐶0 𝜋 0 1 − 𝜋 3 最尤法 ベイズ統計 この関数が最大になる 𝜋 の値を求める この関数の形（×事前分布）がそのまま事後分布の形になる 03 尤度 24

実際の計算 例 5人がクジを3回ずつひいたところ，それぞれ(1,2,2,1,0)回当たった。 このときクジの当たり確率 𝜋 の事後確率は？ 【尤度関数】 𝐿 𝜋 𝒀, 𝑛 = 3𝐶1 𝜋1 1 − 𝜋 2 × 3𝐶2 𝜋 2 1 − 𝜋 1 × 3𝐶2 𝜋 2 1 − 𝜋 1 × 3𝐶1 𝜋 1 1 − 𝜋 2 × 3𝐶0 𝜋 0 1 − 𝜋 3 ▌尤度関数の積 × 𝐿 𝜋 𝑌1 , 𝑛 = 3𝐶1 𝜋1 1 − 𝜋 2 ＝ 𝐿 𝜋 𝑌2 , 𝑛 = 3𝐶2 𝜋 2 1 − 𝜋 1 03 尤度 × とりあえず最初の２つだけ掛け算してみると… 𝐿 𝜋 𝑌1 , 𝑛 × 𝐿 𝜋 𝑌2 , 𝑛 25

実際の計算 例 5人がクジを3回ずつひいたところ，それぞれ(1,2,2,1,0)回当たった。 このときクジの当たり確率 𝜋 の事後確率は？ 【尤度関数】 𝐿 𝜋 𝒀, 𝑛 = 3𝐶1 𝜋1 1 − 𝜋 2 × 3𝐶2 𝜋 2 1 − 𝜋 1 × 3𝐶2 𝜋 2 1 − 𝜋 1 × 3𝐶1 𝜋 1 1 − 𝜋 2 × 3𝐶0 𝜋 0 1 − 𝜋 3 ▌尤度関数の積 × 𝐿 𝜋 𝑌1 , 𝑛 = 3𝐶1 𝜋1 1 − 𝜋 2 ＝ 𝐿 𝜋 𝑌2 , 𝑛 = 3𝐶2 𝜋 2 1 − 𝜋 1 03 尤度 × ３つの設定（𝜋 = 0.2, 0.5, 0.8）について具体的に計算すると 𝐿 𝜋 𝑌1 , 𝑛 × 𝐿 𝜋 𝑌2 , 𝑛 26

実際の計算 例 5人がクジを3回ずつひいたところ，それぞれ(1,2,2,1,0)回当たった。 このときクジの当たり確率 𝜋 の事後確率は？ 【尤度関数】 𝐿 𝜋 𝒀, 𝑛 = 3𝐶1 𝜋1 1 − 𝜋 2 × 3𝐶2 𝜋 2 1 − 𝜋 1 × 3𝐶2 𝜋 2 1 − 𝜋 1 × 3𝐶1 𝜋 1 1 − 𝜋 2 × 3𝐶0 𝜋 0 1 − 𝜋 3 ▌尤度関数の積 小さくなりすぎるので スケールを変えて表示 同じ要領で５つすべて掛け算してみると… × 1 2 3𝐶1 𝜋 1 − 𝜋 × 3𝐶2 𝜋 2 1−𝜋 1 × 3𝐶2 𝜋 2 1−𝜋 1 ＝ × 1 2 3𝐶1 𝜋 1 − 𝜋 0 3 3𝐶0 𝜋 1 − 𝜋 03 尤度 𝐿 𝜋 𝒀, 𝑛 27

実際の計算 例 5人がクジを3回ずつひいたところ，それぞれ(1,2,2,1,0)回当たった。 このときクジの当たり確率 𝜋 の事後確率は？ 【尤度関数】 𝐿 𝜋 𝒀, 𝑛 = 3𝐶1 𝜋1 1 − 𝜋 2 × 3𝐶2 𝜋 2 1 − 𝜋 1 × 3𝐶2 𝜋 2 1 − 𝜋 1 × 3𝐶1 𝜋 1 1 − 𝜋 2 × 3𝐶0 𝜋 0 1 − 𝜋 3 ▌この中で設定を3通りだけにして考えてみると 設定2 設定1 設定3 03 尤度 28

実際の計算 例 5人がクジを3回ずつひいたところ，それぞれ(1,2,2,1,0)回当たった。 このときクジの当たり確率 𝜋 の事後確率は？ 【尤度関数】 𝐿 𝜋 𝒀, 𝑛 = 3𝐶1 𝜋1 1 − 𝜋 2 × 3𝐶2 𝜋 2 1 − 𝜋 1 × 3𝐶2 𝜋 2 1 − 𝜋 1 × 3𝐶1 𝜋 1 1 − 𝜋 2 × 3𝐶0 𝜋 0 1 − 𝜋 3 ▌この中で設定を3通りだけにして考えてみると 設定1 𝜋 = 0.2 𝐿 𝜋 = 0.2 𝒀, 𝑛 = 0.384 × 0.096 × 0.096 × 0.384 × 0.512 ≃ 0.000696 設定2 𝜋 = 0.5 𝐿 𝜋 = 0.5 𝒀, 𝑛 = 0.375 × 0.375 × 0.375 × 0.375 × 0.125 ≃ 0.002472 設定3 𝜋 = 0.8 𝐿 𝜋 = 0.8 𝒀, 𝑛 = 0.096 × 0.384 × 0.384 × 0.096 × 0.008 ≃ 0.000011 ◀ 最尤法 ならこれ 全データの積を計算するのは大変＆値がかなり小さくなる ▶ 実際には対数に変換したうえで計算を行います 03 尤度 29

対数尤度関数 ▌尤度関数 𝑁 はデータの総数 𝑁 𝐿 𝜋 𝒀, 𝑛 = ෑ 𝑃 𝑌𝑖 𝜋, 𝑛 = 𝑃 𝑌1 𝜋, 𝑛 × ⋯ × 𝑃 𝑌𝑁 𝜋, 𝑛 𝑖=1 ▌対数尤度関数 (Log-likelihood function) 𝑁 𝐿𝐿 𝜋 𝒀, 𝑛 = ෍ log 𝑃 𝑌𝑖 𝜋, 𝑛 = log 𝑃 𝑌1 𝜋, 𝑛 + ⋯ + log 𝑃 𝑌𝑁 𝜋, 𝑛 𝑖=1 掛け算から足し算に変化するので計算しやすい！ 03 尤度 30

勝手に対数をとってもいいのか？ ▌p. 26の尤度関数を対数に変換してみると… 尤度関数 対数尤度関数 形は変わるが頂点の場所は変わらない 03 尤度 最尤法 ベイズ統計 このまま使える 事後分布はexp 𝐿𝐿 𝜋 𝑌 = 𝑘 31

パラメータが複数あるとき 例 ある母集団から5人をランダムサンプリングし，身長を測ったところ， 𝑌 = (171, 168, 170, 175, 180) であった。 このとき母集団の身長の平均値 𝜇 および標準偏差 𝜎 はそれぞれいくつ？ 【尤度関数】 𝑁 171 − 𝜇 2 𝐿 𝜇, 𝜎 𝒀 = ෑ 𝑃 𝑌𝑖 𝜇, 𝜎 = exp − 2 2𝜎 2 2𝜋𝜎 𝑖=1 ２変数関数 1 168 − 𝜇 2 1 170 − 𝜇 2 × exp − × exp − 2 2 2 2𝜎 2𝜎 2 2𝜋𝜎 2𝜋𝜎 1 175 − 𝜇 2 1 180 − 𝜇 2 × exp − × exp − 2 2 2 2𝜎 2𝜎 2 2𝜋𝜎 2𝜋𝜎 1 03 尤度 32

あるいは２つまとめて尤度関数 ▌2つのパラメータを同時に動かしながら𝐿(𝜇, 𝜎|𝑌)を計算していくと… 等高線プロット 最尤法 𝑌 = (171, 168, 170, 175, 180) というデータに対して は 𝜇, 𝜎 = 173,4.3 くらいのときが尤度が最大 ▶ 最尤推定値は 𝜇, 𝜎 = 173,4.3 ベイズ統計 これに事前分布をかけたものがそのまま 同時事後確率分布となる 𝑁 𝐿 𝜇, 𝜎 𝒀 = ෑ 𝑃 𝑌𝑖 𝜇, 𝜎 𝑖=1 03 尤度 33

２ stanで尤度関数をつくる 03 尤度 34

ベイズ推定をやってみよう ベイズの定理 ▌ベイズ推定をするためには 当は尤度と事前分布を両方指定する必要がある 𝑃 𝜃𝑌 ∝𝑃 𝑌𝜃 𝑃 𝜃 ▌とりあえず今回は 事前分布を一様分布として（つまり尤度関数だけで）ベイズ推定をやってみます stanの 𝑃 𝜃𝑌 文法を少しずつ説明していきたいと思います ∝ × 𝑃 𝑌𝜃 𝑃(𝜃) 最もシンプルな例として，このケースを試してみましょう。 例 クジを3回ひいたらちょうど1回当たった。このときクジの当たり確率 𝜋 の事後分布は？ 03 尤度 35

MCMCによるベイズ推定 ▌stanによる推定の流れ ① stanコードを書く ② Rでstanコードをコンパイルする（cmdstan_model関数） ③ Rでデータを準備する ④ 実際にサンプリングを行う（$sample関数） ⑤ 結果を見る（いろいろな関数） 本気で覚えたいと思う人はただ書き写すのではなくコードの意味を考えながら書いてください。 大して 気じゃない人はとりあえず写してください。 03 尤度 36

stanコードの基本ルール 拡張子は” stan”にしておくのが通例 実際にはただのテキストファイルなので，メモ帳などで開ける Rstudioではサジェストなどしてくれるので，”.stan”にしておくとやりやすい 各行の最後にはセミコロンを付ける “//”以後はコメント扱い 変数の型を指定する必要がある C++に変換されるので，例えば整数型ならintとか指定しないと動かない 最後の行は空白でおわる 実際には空白じゃなくても動くが，ちょっと注意される 03 尤度 37

Rstudioの場合 ▌左上からファイルを新規作成 色々書いてあるテンプレートが出てくる コメントアウトされているものや 各ブロックの中身を一旦空にしておく 03 尤度 38

stanコードのブロック data { どんな形のデータ(𝑌)が与えられるかを指定する。 ※stanコードでは「こんな形式のデータが来る」を指定します。 実際のデータはRから渡します。 } parameters { 推定するパラメータ(𝜃)を指定する。 } model { 実際に事後分布の形を規定するもの(𝑃 𝑌 𝜃 𝑃 𝜃 )を指定する。 すなわち事前分布𝑃 𝜃 と尤度𝑃 𝑌 𝜃 の両方を書くブロック。 } 𝑃 𝜃𝑌 = 𝑃 𝑌𝜃 𝑃 𝜃 𝑃 𝑌 03 尤度 39

stanコードのブロック data { どんな形のデータ(𝑌)が与えられるかを指定する。今回の例では • 試行数 • 当たり数 ※stanコードでは「こんな形式のデータが来る」を指定します。 実際のデータはRから渡します。 } の2つが与えられている。 parameters { 推定するパラメータ(𝜃)を指定する。 } model { 実際に事後分布の形を規定するもの(𝑃 𝑌 𝜃 𝑃 𝜃 )を指定する。 すなわち事前分布𝑃 𝜃 と尤度𝑃 𝑌 𝜃 の両方を書くブロック。 } 𝑃 𝜃𝑌 = 𝑃 𝑌𝜃 𝑃 𝜃 𝑃 𝑌 03 尤度 40

dataブロック ▌stanとRの関係 ① データ (𝑌) を与える ②𝑃 𝜃 𝑌 ∝𝑃 𝑌 𝜃 𝑃 𝜃 に基づいて 事後分布𝑃 𝜃 𝑌 を生成する ③ 事後分布 (𝑃(𝜃|𝑌)) を返す ▌dataブロックで「どんなデータが与えられるか」を記す 今回は「試行数」と「当たり回数」 成功数はKだろうとXだろうとわかれば良いです（SUCCESSとかでもOK） ただ、入る値は必ず整数である、ということは明示します ※stanの二項分布は「整数 (int) しか受け付けない」ようになっています data { int N; 整数（integer）型であることを明示 int K; } 03 尤度 (変数型) (変数名); が基本形 41

stanコードのブロック data { どんな形のデータ(𝑌)が与えられるかを指定する。今回の例では int N; • 試行数 int K; • 当たり数 ※stanコードでは「こんな形式のデータが来る」を指定します。 実際のデータはRから渡します。 } の2つが与えられている。 parameters { 推定するパラメータ(𝜃)を指定する。今回の例では • 成功確率 } の1つだけ。 model { 実際に事後分布の形を規定するもの(𝑃 𝑌 𝜃 𝑃 𝜃 )を指定する。 すなわち事前分布𝑃 𝜃 と尤度𝑃 𝑌 𝜃 の両方を書くブロック。 } 𝑃 𝜃𝑌 = 𝑃 𝑌𝜃 𝑃 𝜃 𝑃 𝑌 03 尤度 42

parametersブロック ▌今回推定したいのは「当たり確率」 これもthetaだろうとpiだろうとわかれば良いです（もちろん重複は無し） ただ，あり得る値は0から1のみ，ということは明示します ※stanの二項分布の成功率パラメータは「0から1しかとらない」という前提で作られています また、単純に推定を安定させる働きもあります parameters { 実数（real）型であることを明示 数字の型はざっくり分けると • int （整数のみ） • real （実数：小数もあり） • complex （複素数） の３種類です real <lower=0,upper=1> pi; } 取りうる値の範囲を指定する書き方 lowerだけ、upperだけでも可能 一応dataブロックなどでも同じようにかける 03 尤度 43

stanコードのブロック data { どんな形のデータ(𝑌)が与えられるかを指定する。今回の例では int N; • 試行数 int K; • 当たり数 ※stanコードでは「こんな形式のデータが来る」を指定します。 実際のデータはRから渡します。 } の2つが与えられている。 parameters { 推定するパラメータ(𝜃)を指定する。今回の例では real <lower=0,upper=1> pi; • 成功確率 } の1つだけ。 model { 実際に事後分布の形を規定するもの(𝑃 𝑌 𝜃 𝑃 𝜃 )を指定する。 すなわち事前分布𝑃 𝜃 と尤度𝑃 𝑌 𝜃 の両方を書くブロック。 } 𝑃 𝜃𝑌 = 𝑃 𝑌𝜃 𝑃 𝜃 𝑃 𝑌 03 尤度 が，今日だけは尤度のみで 44

modelブロック ▌今回の尤度𝑃 𝑌 𝜃 は二項分布 事前分布は指定しなければ 勝手に一様分布になります 大抵の確率分布はstanに実装されているので、マニュアルを見る ※Rとは違う名前になってることが多いので注意 分布の名称・パラメータの順番は model { “~”は「分布に従う」という意味 「確率的に発生する」とか 「その分布から生じる」とか K ~ binomial(N, pi); } マニュアルを見てください ※ 気でstanを使えるようになりたかっ たら、自分でマニュアルを読めるように なるとよい 一方で代入する場合は”=”を使う （そのうち出てくると思います） 03 尤度 45

完成したstanコード data { どんな形のデータ(𝑌)が与えられるかを指定する。今回の例では int N; • 試行数 int K; • 当たり数 ※stanコードでは「こんな形式のデータが来る」を指定します。 実際のデータはRから渡します。 } の2つが与えられている。 parameters { 推定するパラメータ(𝜃)を指定する。今回の例では real <lower=0,upper=1> pi; • 成功確率 } の1つだけ。 model { 実際に事後分布の形を規定するもの(𝑃 𝑌 𝜃 𝑃 𝜃 )を指定する。 K ~ binomial(N, pi); } すなわち事前分布𝑃 𝜃 と尤度𝑃 𝑌 𝜃 の両方を書くブロック。 𝑃 𝜃𝑌 = 𝑃 𝑌𝜃 𝑃 𝜃 𝑃 𝑌 03 尤度 が，今日だけは尤度のみで 46

登場人物の対応関係 【今回のベイズモデリング】 data { 𝑛 回クジをひいたときの当たり回数 𝑘 は，確率 𝜋 の二項分布に従う 𝐾 ∼ 𝐵𝑖𝑛𝑜𝑚𝑖𝑎𝑙(𝑁, 𝜋) int N; int K; } parameters { real <lower=0,upper=1> pi; } model { K ~ binomial(N, pi); } 【登場人物の紹介】 dataブロックとparametersブロックの中で モデリングの登場人物を紹介する …今回は𝐾, 𝑁, 𝜋の３つだけ ▶データ等によって所与のものはdataブロックに書き， データじゃない（推定したい）ものはparametersブロックに書く 【modelブロック】 登場人物の関係性を書いてあげる …と言っても上のベイズモデリングをstanの記法で書き写すだけ ▶ここでdataブロックとparametersブロックのいずれにも 登場していないものが出てくるとおかしい 03 尤度 47

plate notation （プレート記法） 【今回のベイズモデリング】 data { 𝑛 回クジをひいたときの当たり回数 𝑘 は，確率 𝜋 の二項分布に従う 𝐾 ∼ 𝐵𝑖𝑛𝑜𝑚𝑖𝑎𝑙(𝑁, 𝜋) int N; int K; ▌ plate notationとは } データとパラメータの関係を視覚的に表現する方法の一つ parameters { real <lower=0,upper=1> pi; 【最低限の決まりごと】 (Lee & Wagenmakers, 2013 による) • 潜在変数（パラメータ）は白い丸 • 観測変数（データなど）はグレーの四角 } model { で表す 【今回のベイズモデリングのplate notation】 𝑁 K ~ binomial(N, pi); } 𝐾 𝜋 03 尤度 48

MCMCによるベイズ推定 ▌stanによる推定の流れ ① stanコードを書く 出来上がったらstanファイルを保存します この資料では ”model_binom.stan” という名前で保存したとして進めます ② Rでstanコードをコンパイルする（cmdstan_model関数） ③ Rでデータを準備する ④ 実際にサンプリングを行う（$sampling関数） ⑤ 結果を見る（いろいろな関数） 03 尤度 49

MCMCによるベイズ推定 ▌stanによる推定の流れ ① stanコードを書く ② Rでstanコードをコンパイルする（cmdstan_model関数） library(cmdstanr) # Rを開くたびに毎回行ってください model <- cmdstan_model(“model_binom.stan”) 細かい説明は省略しますが 適当な名前のオブジェクトに代入する必要あり ③ Rでデータを準備する ④ 実際にサンプリングを行う（$sample関数） ⑤ 結果を見る（いろいろな関数） 03 尤度 50

MCMCによるベイズ推定 ▌stanによる推定の流れ ① stanコードを書く ② Rでstanコードをコンパイルする（cmdstan_model関数） ③ Rでデータを準備する stan_data <- list(N=3, K=1) list型は「どんな形のものでもとりあえず 【ポイント】 くっつけて一つのオブジェクトにできる」型 • list型で用意する • stanコードのdataブロックで指定したものをすべて含む ④ 実際にサンプリングを行う（$sample関数） ⑤ 結果を見る（いろいろな関数） 03 尤度 51

MCMCによるベイズ推定 ▌stanによる推定の流れ ① stanコードを書く ② Rでstanコードをコンパイルする（cmdstan_model関数） ③ Rでデータを準備する ④ 実際にサンプリングを行う（$sample関数） result <- model$sample(data = stan_data) 他にもいろいろなオプションはあるのですが それについてはおいおい説明します… All 4 chains finished successfully. Mean chain execution time: 0.0 seconds. Total execution time: 0.7 seconds. 的な表示が出たら無事成功！ ⑤ 結果を見る（いろいろな関数） 03 尤度 52

MCMCによるベイズ推定 ▌stanによる推定の流れ ① stanコードを書く ② Rでstanコードをコンパイルする（cmdstan_model関数） ③ Rでデータを準備する ④ 実際にサンプリングを行う（$sample関数） ⑤ 結果を見る（いろいろな関数） いろいろな方法があるのですが，とりあえず今回は事後分布の形を見てみましょう 03 尤度 53

事後分布の形を見る ▌追加パッケージの準備 cmdstanrはあくまでもMCMCによる乱数生成を行うだけ その結果をもとにプロットしたりする場合にはもうひと手間必要 これを楽にやってくれるパッケージの一つがbayesplotパッケージ install.packages(“bayesplot”) library(bayesplot) ▌早速事後分布を書いてみる mcmc_dens(result$draws(), pars = "pi") bayesplotパッケージには mcmc_***()という名前の関数が色々あります （他の関数はおいおい説明します…） result$draws()をすると， MCMCで生成された乱数を取り出すことが出来ます 03 尤度 54

事後分布はどんな形？ 【解析的に求めた尤度関数】 mcmc_dens(result$draws(), pars = "pi") 𝑃 𝜃 = 𝜋 𝑘 = 1 = 3𝐶1 𝜋 1 1 − 𝜋 2 乱数なので微妙に結果が変わるかもしれませんが 概ね同じような形になっているはずです 確かに大体同じ形を 推定することが出来ました 03 尤度 55

まとめと次回予告 【まとめ】 ▌ベイズ統計学における尤度関数の意味が分かりました 的には尤度関数がそのまま事後分布を形作っている ▌stanの第一歩を踏み出しました 細かいことはともかく， 的な流れを説明しました 【次回予告】 そろそろ事前分布のことを考えてあげましょう stanで推定できた事後分布から点推定・区間推定などしてみましょう 03 尤度 56