ベイズ統計_03_尤度

5.7K Views

April 23, 24

スライド概要

神戸大学大学院経営学研究科で2024年度より開講している「ベイズ統計」の講義資料「03_尤度」です。尤度および尤度関数についておさらいした後で,ベイズ統計学における尤度関数の使い方を説明しています。また,stanで尤度関数がうまく復元できることを確認しつつstanの基本中の基本的なことを説明しています。

profile-image

神戸大学経営学研究科准教授 分寺杏介(ぶんじ・きょうすけ)です。 主に心理学的な測定・教育測定に関する研究を行っています。 講義資料や学会発表のスライドを公開していきます。 ※スライドに誤りを見つけた方は,炎上させずにこっそりお伝えいただけると幸いです。

シェア

またはPlayer版

埋め込む »CMSなどでJSが使えない場合

関連スライド

各ページのテキスト
1.

ベイズ統計 03 尤度 分寺 杏介 神戸大学大学院 経営学研究科  [email protected] ※本スライドは,クリエイティブ・コモンズ 表示-非営利 4.0 国際 ライセンス(CC BY-NC 4.0)に従って利用が可能です。

2.

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

3.

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

4.

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

5.

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

6.

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

7.

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

8.

1 尤度とは 03 尤度 8

9.

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

10.

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

11.

各設定のもとで「データが得られる確率」を計算 • 以下のうち,今日の 設定 はどれでしょうか? 設定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

12.

連続変数に拡張してみると 例 クジを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

13.

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

14.

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

15.

確率(密度)関数を別の見方で パラメータが所与の場合(確率関数) ▌𝑛 = 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

16.

ただし「確率」ではない パラメータが所与の場合(確率関数) ▌𝑛 = 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

17.

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

18.

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

19.

前回登場した面積による考え方 ▌「すべてのあり得る 設定 から,どれが最もありえそうか」を判断 設定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

20.

前回登場した面積による考え方 ▌「すべてのあり得る 設定 から,どれが最もありえそうか」を判断 設定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

21.

事前の信念による調整もあり ▌(例)過去の傾向から,今日は 設定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

22.

事後分布は「確率分布」なので… ベイズの定理 𝑃 𝑌𝜃 𝑃 𝜃 𝑃 𝜃𝑌 = 𝑃 𝑌 【実際に 設定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

23.

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

24.

複数のデータがあるときの尤度関数 例 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

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 尤度 × とりあえず最初の2つだけ掛け算してみると… 𝐿 𝜋 𝑌1 , 𝑛 × 𝐿 𝜋 𝑌2 , 𝑛 25

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 , 𝑛 = 3𝐶1 𝜋1 1 − 𝜋 2 = 𝐿 𝜋 𝑌2 , 𝑛 = 3𝐶2 𝜋 2 1 − 𝜋 1 03 尤度 × 3つの設定(𝜋 = 0.2, 0.5, 0.8)について具体的に計算すると 𝐿 𝜋 𝑌1 , 𝑛 × 𝐿 𝜋 𝑌2 , 𝑛 26

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 ▌尤度関数の積 小さくなりすぎるので スケールを変えて表示 同じ要領で5つすべて掛け算してみると… × 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

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通りだけにして考えてみると 設定2 設定1 設定3 03 尤度 28

29.

実際の計算 例 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

30.

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

31.

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

32.

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

33.

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

34.

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

35.

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

36.

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

37.

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

38.

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

39.

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

40.

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

41.

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

42.

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

43.

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

44.

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

45.

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

46.

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

47.

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

48.

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

49.

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

50.

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

51.

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

52.

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

53.

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

54.

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

55.

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

56.

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