ベイズ統計_05_基本的なベイズ推論(2)
13.2K Views
May 21, 24
スライド概要
神戸大学大学院経営学研究科で2024年度より開講している「ベイズ統計」の講義資料「05_基本的なベイズ推論2」です。前の資料に引き続き，基本的な問題としてポアソン分布・正規分布のパラメータ推定問題について，最尤法・共役事前分布を用いた解析的な導出・stanでの推定を行っていきます。
神戸大学経営学研究科准教授 分寺杏介（ぶんじ・きょうすけ）です。 主に心理学的な測定・教育測定に関する研究を行っています。 講義資料や学会発表のスライドを公開していきます。 ※スライドに誤りを見つけた方は，炎上させずにこっそりお伝えいただけると幸いです。
関連スライド
各ページのテキスト
ベイズ統計 05 基本的なベイズ推論(2) 分寺 杏介 神戸大学大学院 経営学研究科 [email protected] ※本スライドは，クリエイティブ・コモンズ 表示-非営利 4.0 国際 ライセンス（CC BY-NC 4.0）に従って利用が可能です。
（前回のおさらい）事前分布の決め方 1. そもそもなぜパラメータ𝜋の事前分布にベータ分布を置いたの？ 共役事前分布だから（ただしstanの中では別に一様分布とかでもOK） 2. じゃあ他の場合はどんな分布を置けばいいの？ ケースバイケースなので，先行研究を参考にしたりしましょう。 今後の授業でも少しずつ分かってくると思います。 3. 事前分布の形はどうやって決めたら怒られないの？ 05 基本的なベイズ推論(2) 2
（前回のおさらい）Gelmanのおすすめ（抜粋） Gelmanって誰？ コロンビア大学の統計学の教授 stanを作った人 Bayesian Data Analysis (BDA) の第一著者 事前分布の５つのレベル 標準化されたスケールでの話 基本的に事前分布は「どの程度の情報を含むか」がカギ • Flat prior (not usually recommended); • Super-vague but proper prior: normal(0, 1e6) (not usually recommended); • Weakly informative prior, very weak: normal(0, 10); • Generic weakly informative prior: normal(0, 1); • Specific informative prior: normal(0.4, 0.2) or whatever. 05 基本的なベイズ推論(2) 3
（前回のおさらい）Gelmanのおすすめ（抜粋） 弱情報事前分布の場合は どう考えてもありえない値は出ないように もしかしたらありえる値は出てもおかしくないように Fully informativeよりはWeakly informativeに 事前分布の情報が減ると基本的には推定精度が下がるが， それよりも広い範囲の値をカバーできる頑健性のメリットのほうが大きい When using informative priors, be explicit about every choice ベイズ統計では尤度と事前分布の両方 を分析者が決める必要があるので 事前分布を報告しないと 一様分布や切断分布よりは弱情報事前分布がよい と思われてしまいます！ 常に使用した事前分布は明示すること 「コイツ分かってないな」 値域が明確に決まる場合でないと制約がきつすぎる 05 基本的なベイズ推論(2) 4
ということで今日も 実際の分析をやっていきたいと思います 基本的には同じ分析を ① 非ベイズ的方法（最尤法や標本理論的な仮説検定） 時間がなさそうだったらスキップ ② 解析的なベイズ（共役事前分布を利用した方法） ③ 数値計算的なベイズ（StanでMCMC） の1+2パターンでやっていきます。 その過程で stanに慣れましょう できればRにも慣れましょう 05 基本的なベイズ推論(2) 5
１ 分析実践編(2) ポアソン分布のパラメータ 05 基本的なベイズ推論(2) 6
事例 例 大地震は怖いので、防災対策をしようと思いました。 でも、そもそも地震がすぐには起こらないなら対策の必要もありません。あなたは各年に発生 した大地震の件数から、１年当たり平均何件大地震が起こるかを推定することにしました。 これを使えば10年以内に大地震が起こる確率を予測できます。 【データの読み込み】 この名前はなんでもOK "data_quake.csv"をワーキングディレクトリに配置して dat <- read.csv("data_quake.csv") 中身はこんなデータ year:何年のデータか (1919から2023年までの各年 105行) time:その年に何回の大地震があったか (震度5弱以上の回数) ※データは気象庁から取得可能です 05 基本的なベイズ推論(2) 7
事例 例 大地震は怖いので、防災対策をしようと思いました。 でも、そもそも地震がすぐには起こらないなら対策の必要もありません。あなたは各年に発生 した大地震の件数から、１年当たり平均何件大地震が起こるかを推定することにしました。 これを使えば10年以内に大地震が起こる確率を予測できます。 【まずはなにより事例の整理】 推定に必要な情報 今回の事例 データ 𝑌 各年におきた地震の数 𝑥 = (0,1,1,3,21,3, ⋯ ) 推定したいパラメータ 𝜃 年間の平均地震発生件数 𝜆 尤度 𝑃(𝑌|𝜃) ポアソン分布 𝑃𝑜𝑖𝑠(𝑥|𝜆) 事前分布 𝑃(𝜃) ガンマ分布 𝐺𝑎𝑚𝑚𝑎(𝜆|𝛼, 𝛽) こいつが今回の共役事前分布 A パラメータの事後分布は ガンマ 分布になる 05 基本的なベイズ推論(2) 8
カウントデータに対する代表的な確率分布 シメオン・ドニ・ポアソン(1781-1840) Wikipediaより ポアソン過程に基づく事象の回数 確率的に発生する事象が一定時間のうちに何回起こったか ただしその事象は基本的に起こりにくい 各事象はほかの事象とは無関係 「一定時間」のフレームを小さくした場合ほぼ０ということ ある事象が起こりやすくなる確変的なものや時系列的なことは考えない ポアソン分布に合う事象の例 １日にコンビニに来る客の数 １ヶ月に起こる地震の数 １年に倒産する企業の数 歴史的には「馬に蹴られて死亡する兵士の数」 といったものにもポアソン分布がよく当てはまりました 05 基本的なベイズ推論(2) 9
ポアソン分布｜Poisson distribution 𝜆𝑥 𝑒 −𝜆 𝑃 𝑋 = 𝑥|𝜆 = 𝑥! 関数 パラメータ 𝜆 𝑒 = 2.718 … （ネイピア数） 𝜋 みたいな特殊な定数 一定時間に起こる回数の期待値 略記 𝑃𝑜𝑖𝑠(𝜆) 平均値 𝜆 分散 𝜆 平均値と分散が同じ という特徴をもつ パラメータを変えると様々な形の関数が書ける 𝑃𝑜𝑖𝑠(0.5) 𝑃𝑜𝑖𝑠(2) 05 基本的なベイズ推論(2) 𝑃𝑜𝑖𝑠(10) 10
ポアソン分布の計算例 問 Aさんはスマホを半年に2回くらい落としてしまいます。 では，Aさんがスマホを１年に１回も落とさない確率はいくつでしょうか。 𝜆𝑥 𝑒 −𝜆 𝑃 𝑋=𝑥 = 𝑥! ポアソン分布の性質 時間のフレームを変えると 𝜆 が変わる 「Aさんが半年のうちにスマホを落とす回数」の確率分布は𝑃𝑜𝑖𝑠(2) 「Aさんが１年のうちにスマホを落とす回数」の確率分布は𝑃𝑜𝑖𝑠(4)となる 40 𝑒 −4 𝑃𝑜𝑖𝑠(4) に𝑥 = 0 を入れると 𝑃 𝑋 = 0|𝜆 = 4 = = 𝑒 −4 ≃ 0.0183 0! 05 基本的なベイズ推論(2) 11
ガンマ分布｜Gamma distribution 何かしらの出来事が発生するまでの時間に関する確率分布 ある期間中に平均 𝛽 回発生する出来事が 𝛼 回起こるまでの時間 𝛽 𝛼 𝛼−1 −𝛽𝑥 𝑃 𝑥 𝛼, 𝛽 = 𝑥 𝑒 Γ 𝛼 関数 パラメータ 略記 期待値 分散 𝛼 ある出来事の回数 𝛽 ある出来事の発生頻度 Γ 𝛼 は「ガンマ関数」ですが 正規化定数の一部なので気にしなくてOKです 【当てはまる例】 𝐺𝑎𝑚𝑚𝑎(𝑥|𝛼, 𝛽) 𝛼 𝛽 𝛼 𝛽2 05 基本的なベイズ推論(2) 𝛽 スマホを年３回くらい落とす人が ２回落とすまでの時間（年） 𝛼 𝐺𝑎𝑚𝑚𝑎(𝑥|2, 3) 𝛽 10年に1度起きる大地震が １回起きるまでの時間（年） 𝛼 𝐺𝑎𝑚𝑚𝑎(𝑥|1, 0.1) 12
ガンマ分布の特徴 0以上の値を取る 「時間」を表す分布なので当然 パラメータによって多様な形に変わる 正の値をとる変数に対する事前分布として 用いられることが多い 𝛼 = 2のとき 𝛽=5 𝛽=3 𝛽=1 指数分布に従う確率変数の和でもある 指数分布は「ある期間中に平均 𝛽 回発生する 出来事が１回起こるまでの時間」 これを 𝛼 個足したものは 𝐺𝑎𝑚𝑚𝑎(𝑥|𝛼, 𝛽) 𝛽 = 1のとき 𝛼=5 𝛼=3 𝛼=1 つまりガンマ分布では，𝛼 回の事象は独立に生じる （前回の発生にどれだけの時間がかかったかは無関係） であることが仮定されている 05 基本的なベイズ推論(2) 13
（補足）共役事前分布の見つけ方 カーネルの形から考える 事後分布は尤度×事前分布の形に比例するので，今回の場合 𝜆𝑥 𝑒 −𝜆 𝑃 𝜃𝑌 = × 𝑃(𝜃) 𝑥! 尤度＝ポアソン分布 のカーネルは 𝜆𝑥 𝑒 −𝜆 正規化定数 × 𝜆？ 𝑒 ？𝜆 の形ならば 尤度をかけた後でも同じ 正規化定数 × 𝜆？ 𝑒 ？𝜆 の形になる！ 事前分布𝑃(𝜆)が ガンマ分布 𝛽𝛼 𝛼−1 −𝛽𝜆 𝑃 𝜆 𝛼, 𝛽 = 𝜆 𝑒 Γ 𝛼 05 基本的なベイズ推論(2) 14
非ベイズ的点推定（最尤推定） ポアソン分布の尤度＝確率質量関数 𝐿 𝜆𝑥 =𝑓 𝑥𝜆 1919年 1920年 𝜆𝑥1 𝑒 −𝜆 まず1年目については𝐿 𝜆 𝑥1 = ，2年目は𝐿 𝜆 𝑥2 𝑥1 ! 𝜆𝑥 𝑒 −𝜆 = 𝑥! 𝜆𝑥2 𝑒 −𝜆 = ，… 𝑥2 ! この調子で105年分 各年の地震発生件数が独立であると仮定すると 𝑛 𝜆𝑥𝑖 𝑒 −𝜆 𝜆sum(𝒙) × 𝑒 −𝑛𝜆 𝐿 𝜆𝒙 =ෑ = ς𝑛𝑖=1 𝑥𝑖 ! 𝑥𝑖 ! 𝑛 sum 𝒙 = 𝑖=1 𝑖=1 この関数が最大値を取る 𝜆 の値を探します 05 基本的なベイズ推論(2) 𝑥𝑛 15
非ベイズ的点推定 𝑛 𝜆𝑥𝑖 𝑒 −𝜆 𝜆sum 𝒙 × 𝑒 −𝑛𝜆 𝐿 𝜆𝒙 =ෑ = ς𝑛𝑖=1 𝑥𝑖 ! 𝑥𝑖 ! 𝑖=1 対数とって 𝜆 で微分して 𝑛 log𝐿 𝜆 𝒙 = log 𝜆 sum 𝒙 − 𝑛𝜆 − log ෑ 1 log𝐿 𝜆 𝒙 = sum 𝒙 − 𝑛 𝜆 これがゼロになる点は ′ 【最尤推定値】 1 𝜆 = sum 𝒙 𝑛 05 基本的なベイズ推論(2) 𝑥𝑖 ! 𝑖=1 カーネルは尤度関数の形と無関係（𝜆が無い） なので，必ず計算途中で消える つまり標本平均ということ 16
非ベイズ的点推定 ということでデータの平均を求めてみる mean(dat$time) [1] 5.619048 大地震は1年間に平均5.619回起こる！ 1年間に発生する 大地震の回数の確率分布 𝑃𝑜𝑖𝑠(𝑥|𝜆 = 5.619) 【問】 このポアソン分布に従うとき， 日本のどこかで1年のうちに1回でも大地震が発生する確率は？ 【答】 1-dpois(0, 5.619) [1] 0.9963717 およそ99.64% 05 基本的なベイズ推論(2) 17
非ベイズ的区間推定 （いくつか方法がありますが…一つの方法をご紹介します） 中心極限定理を使えば楽ちん ポアソン分布𝑃𝑜𝑖𝑠(𝜆)は平均 𝜆，分散𝜆の確率分布 毎年の発生件数がクジで決まっているとしたら 105回のクジの平均値は確率的に変動するはず 中心極限定理により，𝑛年の（標本）平均は正規分布𝑁 𝜆, 𝜆/𝑛 に近づいていく 105年間の平均発生件数の分布 各年の発生件数の分布 𝑥 ∼ 𝑃𝑜𝑖𝑠(𝑥|𝜆 = 5.619) 𝑥ҧ ∼ 𝑁 5.619, 05 基本的なベイズ推論(2) 5.619 105 18
非ベイズ的区間推定の手順 1 とりあえず95%区間を作る 求めたい区間の上限・下限をそれぞれ𝜆𝐿 , 𝜆𝑈 とする 𝜆𝐿 , 𝜆𝑈 をどのように設定すると 𝑃 𝜆𝐿 ≤ 𝜆 ≤ 𝜆𝑈 = 0.95 となるかを求めたら良い 𝜆𝐿 , 𝜆𝑈 の値をどのように設定したら「 𝜆𝐿 から 𝜆𝑈 の間に真の平均値 𝜆 が含まれている確率（割合）が95%になる」のかを求めたい 2 既知の確率分布に従う統計量になるように変形する 中心極限定理により，標本平均の標本分布は 𝑥ҧ ∼ 𝑁 𝜆, 𝜆/𝑛 で近似できる これを標準化した 𝑧ҧ = ҧ 𝑥−𝜆 𝜆/𝑛 は，標準正規分布に従う 𝑧ҧ = 𝑥ҧ − 𝜆 𝜆/𝑛 ∼ 𝑁(0,1) ҧ 𝑃 𝜆𝐿 ≤ 𝜆 ≤ 𝜆𝑈 の真ん中が 𝑥−𝜆 になるように変形させると 𝜆/𝑛 𝑃 𝑥ҧ − 𝜆𝑈 𝜆/𝑛 ≤ 𝑥ҧ − 𝜆 𝜆/𝑛 ≤ 𝑥ҧ − 𝜆𝐿 𝜆/𝑛 05 基本的なベイズ推論(2) 一旦逆になりますが気にしない 19
非ベイズ的区間推定の手順（つづき） 3 もう一つ95%区間を作る 𝑧ҧ = ҧ 𝑥−𝜆 𝜆/𝑛 が標準正規分布に従う，ということは 𝑃 −1.96 ≤ 𝑧ҧ = ҧ 𝑥−𝜆 𝜆/𝑛 ≤ 1.96 = 0.95 と分かる 終わりだよ… 4 2つの式を対応させると… 3 より 2 より 𝑃 −1.96 ≤ 𝑃 𝑥ҧ − 𝜆𝑈 𝜆/𝑛 ≤ 𝑥ҧ − 𝜆 𝜆/𝑛 𝑥ҧ − 𝜆 𝜆/𝑛 ≤ 1.96 = 0.95 ≤ 𝑥ҧ − 𝜆𝐿 𝜆/𝑛 = 0.95 このままだと真値𝜆がわからないと 信頼区間が作れない… 𝜆𝐿 = 𝑥ҧ − 1.96 𝜆/𝑛 𝜆𝑈 = 𝑥ҧ + 1.96 𝜆/𝑛 𝑥ҧ − 𝜆𝑈 𝜆/𝑛 = −1.96 𝑥ҧ − 𝜆𝐿 𝜆/𝑛 = 1.96 05 基本的なベイズ推論(2) 20
非ベイズ的区間推定の手順（つづき） 5 𝑛 が十分に大きければ 𝜆 のかわりに 𝜆መ を用いて区間を作ったとしても 推定値で代用する 同じ95%の割合で真値 𝜆 を含む区間になる ということです。 1 መ 標本平均の最尤推定量 𝜆 = sum 𝒙 は一致性を持つ 𝑛 サンプルサイズが大きければ 𝜆መ は母数 𝜆 に一致する！ 𝜆መ で置き換え መ 𝜆𝑈 = 𝑥ҧ − 1.96 𝜆/𝑛 𝜆𝑈 = 𝑥ҧ + 1.96 𝜆/𝑛 𝜆መ =5.619 𝑛 = 105 𝜆መ 𝜆መ 𝑃 𝑥ҧ − 1.96 ≤ 𝜆 ≤ 𝑥ҧ + 1.96 𝑛 𝑛 መ 𝜆𝐿 = 𝑥ҧ − 1.96 𝜆/𝑛 𝜆𝐿 = 𝑥ҧ − 1.96 𝜆/𝑛 5.619 − 1.96 今回のデータを 当てはめると = 0.95 すべての標本でこの区間を作った場合 95%の割合で真値 𝜆 が含まれる 5.619 5.619 ≤ 𝜆 ≤ 5.619 + 1.96 105 105 【答】 およそ5.166から6.072 05 基本的なベイズ推論(2) 21
ベイズ推定（まずは解析的に） ポアソン分布のパラメータ𝜆を手計算で行うには共役事前分布が使える 推定に必要な情報 今回の事例 データ 𝑌 各年におきた地震の数 𝑥 = (0,1,1,3,21,3, ⋯ ) 推定したいパラメータ 𝜃 年間の平均地震発生件数 𝜆 尤度 𝑃(𝑌|𝜃) ポアソン分布 𝑃𝑜𝑖𝑠(𝑥|𝜆) 事前分布 𝑃(𝜃) ガンマ分布 𝐺𝑎𝑚𝑚𝑎(𝜆|𝛼, 𝛽) 各年の地震発生件数が独立であると仮定すると，データ全体の尤度は 𝑛 𝜆𝑥𝑖 𝑒 −𝜆 𝜆sum(𝒙) × 𝑒 −𝑛𝜆 𝐿 𝜆𝒙 =ෑ = ς𝑛𝑖=1 𝑥𝑖 ! 𝑥𝑖 ! 𝑖=1 05 基本的なベイズ推論(2) 22
事後分布を求める 𝑃 𝜃 𝑌 ∝ 𝑃 𝑌 𝜃 𝑃 𝜃 について 𝑃 𝜃 にガンマ分布（𝐺𝑎𝑚𝑚𝑎(𝜆|𝛼, 𝛽)）， 𝑃 𝑌 𝜃 にポアソン分布（𝑃𝑜𝑖𝑠(𝑥|𝜆)）を置く 事前分布 尤度 𝜆sum 𝒙 × 𝑒 −𝑛𝜆 𝛽𝛼 𝛼−1 −𝛽𝜆 𝐿(𝜆|𝒙) × 𝑃 𝜆|𝛼, 𝛽 = × 𝜆 𝑒 𝑛 ς𝑖=1 𝑥𝑖 ! Γ 𝛼 𝛽𝛼 = 𝑛 × 𝜆sum 𝒙 + 𝛼−1 × 𝑒 − 𝛽+𝑛 𝜆 ς𝑖=1 𝑥𝑖 ! Γ 𝛼 = 正規化定数 × 𝜆sum 𝒙 + 𝛼−1 × 𝑒 − 𝛽+𝑛 𝜆 事前分布 事後分布 𝐺𝑎𝑚𝑚𝑎 𝛼, 𝛽 が尤度によって更新され𝐺𝑎𝑚𝑚𝑎 𝛼 + sum 𝒙 , 𝛽 + 𝑛 になった 05 基本的なベイズ推論(2) 23
ガンマ分布の更新とパラメータの解釈 𝐺𝑎𝑚𝑚𝑎 𝛼, 𝛽 が尤度によって更新され𝐺𝑎𝑚𝑚𝑎 𝛼 + sum 𝒙 , 𝛽 + 𝑛 になった 𝐺𝑎𝑚𝑚𝑎 0.1, 0.01 𝐺𝑎𝑚𝑚𝑎 590.1, 105.01 • データ0.01年分くらいの強さの信念 0.01年(事前)＋105年(データ)＝105.01年 • 1年に平均10回くらいだろう 0.1回(事前)＋590回(データ)＝590.1回 例 ＝0.01年に平均0.1回くらい データ 尤度 105年で 計590回 05 基本的なベイズ推論(2) 24
信念の強さの比較 【事前の信念】 予想 【データ】 105年で590回発生した 1年に平均10回くらいだろう ポアソン分布の 𝜆 = 10 と予想 自信 「正直言って全く自信はないです」 順当に行けば 590 ポアソン分布の 𝜆 = と予想 105 105年分のデータがあるので データ0.01個分（ほぼゼロ） データ105個分 事前の信念はほぼ自信なし（無情報に近い事前分布）だったので 推論の結果はほぼ完全にデータ（尤度）によって決定されました。 05 基本的なベイズ推論(2) 25
更新前後の期待値 事前分布 事後分布 𝐺𝑎𝑚𝑚𝑎 𝛼, 𝛽 が尤度によって更新され𝐺𝑎𝑚𝑚𝑎 𝛼 + sum 𝒙 , 𝛽 + 𝑛 になった 𝛼 ガンマ分布𝐺𝑎𝑚𝑚𝑎 𝛼, 𝛽 の期待値は 𝛽 事後分布の期待値は パラメータ 意味 事前 分布 𝛽 𝛼 𝛽 事前情報の総量 尤度 𝑛 データの総量 sum 𝒙 𝑛 データの平均値 𝛼 + sum 𝒙 𝛼 sum 𝒙 = + 𝛽+𝑛 𝛽+𝑛 𝛽+𝑛 𝛽 𝛼 𝑛 sum 𝒙 = + 𝛽 +𝑛𝛽 𝛽 +𝑛 𝑛 合計1 事前期待値 事後期待値 = 事前情報のウェイト × 事前期待値 + データのウェイト × (データの平均値) 05 基本的なベイズ推論(2) 26
事前の信念の強さを変えると… 𝐺𝑎𝑚𝑚𝑎 𝛼, 𝛽 が尤度によって更新され𝐺𝑎𝑚𝑚𝑎 𝛼 + sum 𝒙 , 𝛽 + 𝑛 になった 例 𝐺𝑎𝑚𝑚𝑎 1000, 100 𝐺𝑎𝑚𝑚𝑎 1590, 205 • データ100年分くらいの強さの信念 100年(事前)＋105年(データ)＝205年 • 1年に平均10回くらいだろう 1000回(事前)＋590回(データ)＝1590回 ＝100年に平均1000回くらい データ 尤度 データの平均 事前平均 105年で 計590回 05 基本的なベイズ推論(2) 27
信念の強さの比較 【事前の信念】 予想 【データ】 105年で590回発生した １年に平均10回くらいだろう ポアソン分布の𝜆 = 10と予想 自信 「かなり自信あるんです」 「江戸時代の古文書から 100年分のデータ見つけました」 順当に行けば 590 ポアソン分布の 𝜆 = と予想 105 105年分のデータがあるので データ105個分 データ100個分 事前の信念の自信がデータの個数とほぼ同じだったので 推論の結果は両者の予想のほぼ中間になりました。 05 基本的なベイズ推論(2) 28
事前分布のパラメータ設定 ガンマ分布の場合，第2パラメータ 𝛽 が信念の強さを表している 𝛽 の値を小さくするほど無情報事前分布に近づいていく 𝛼 ガンマ分布の分散は 2 𝛽 𝛼 の値は結構重要かもしれない 𝛼 ≤ 1 の場合 （𝛼 = 0.001, 𝛽 = 0.001） 𝛼 > 1 の場合 （𝛼 = 2, 𝛽 = 0.001） 0のところが頂点になる それ以降はほぼフラット 頂点は0以外の場所になる かなり広くて薄い山の形 05 基本的なベイズ推論(2) 29
事後分布の点推定値 これがパラメータの分布だとすると，この分布の代表値を使えば良さそう 主な代表値は３種類 𝐺𝑎𝑚𝑚𝑎 590.1, 105.01 ① 平均値（期待値） 事後期待値（EAP） ② 中央値 事後中央値（MED） ③ 最頻値 事後確率最大推定値（MAP） 05 基本的なベイズ推論(2) 30
点推定値①平均値 𝛼 𝛽 ガンマ分布の期待値は ① 平均値（期待値） 590.1 105.01 EAPは 事後期待値（EAP: Expected A Posteriori） 期待値なので正確には すべてのありえる𝜃について 事後確率𝑃 𝜃 𝑌 による 重み付け平均を取る Rの組み込み関数でやる場合 integrate(function(l){dgamma(l,590.1,105.01)*l},lower=0,upper=10) න 𝜃𝑃 𝜃 𝑌 𝑑𝜃 Θ 今回の場合 整理すると න 𝜆590.1 𝑒 −105.01𝜆 𝑑𝜆 𝜆 5.619 もうちょっと簡単には 1. 右の分布から乱数を作りまくる 2. 平均をとる mean(rgamma(100000,590.1, 105.01)) 05 基本的なベイズ推論(2) 31
点推定値②中央値 ② 中央値 事後中央値（MED: median） 中央値なので正確には qgamma(0.5, 590.1, 105.01) もうちょっと簡単に 事後分布が解析的に わからない場合 5.615 1. 右の分布から乱数を作りまくる 2. 中央値をとる median(rgamma(100000, 590.1, 105.01)) 05 基本的なベイズ推論(2) 32
点推定値③最頻値 ③ 最頻値 事後確率最大推定値（MAP: Maximum A Posteriori） つまりこれは最尤法だ！ Rの組み込み関数でやる場合 optimize(¥(x) dgamma(x,590.1,105.01), interval = c(0,10), maximum = TRUE) もうちょっと簡単に？ 先程までと同じように乱数の最頻値を取ってもだめ（連続変数の場合） 05 基本的なベイズ推論(2) 5.609 33
事後分布からの区間推定 事後分布において，頻度主義的な信頼区間と同じように「範囲」を考えると？ ① 左右の端を切る Equal-tailed interval ② 上から数えていく 𝐺𝑎𝑚𝑚𝑎 590.1, 105.01 Highest posterior density interval ベイズ統計における区間 確信区間または信用区間 と呼ばれます credible interval 05 基本的なベイズ推論(2) 34
ベイズ的区間推定①左右の端を切る ① 左右の端を切る 等裾事後確信区間 (ETI: Equal-tailed interval) 非ベイズ的信頼区間と同じように左右の端を2.5%ずつ切り取ると95% 𝐺𝑎𝑚𝑚𝑎 590.1, 105.01 右の場合の95%確信区間は 5.175 qgamma(0.025,590.1,105.01) から 6.082 の間 qgamma(0.975,590.1,105.01) 05 基本的なベイズ推論(2) 35
ベイズ的区間推定②上から数えていく ② 上から数えていく 最大事後密度確信区間 (HDI: Highest posterior density interval) 確率密度的に「最もありえそうな上位95%」を集めていけば… ありえそうランキング上位95% ▼ 右の場合の95%HDIはだいたい 5.161 から 6.072 の間 05 基本的なベイズ推論(2) 36
stanコードを書いていこう A data { どんな形のデータ(𝑌)が与えられるかを指定する。今回の例では • 各年に起きた大地震の数（105年分） が与えられている。 ※stanコードでは「こんな形式のデータが来る」を指定します。 実際のデータはRから渡します。 } A parameters { 推定するパラメータ(𝜃)を指定する。 } A model { 実際に事後分布の形を規定するもの(𝑃 𝑌 𝜃 𝑃 𝜃 )を指定する。 すなわち事前分布𝑃 𝜃 と尤度𝑃 𝑌 𝜃 の両方を書くブロック。 } 𝑃 𝜃𝑌 = 𝑃 𝑌𝜃 𝑃 𝜃 𝑃 𝑌 05 基本的なベイズ推論(2) 37
dataブロック 今回のデータは「各年の地震の数」 以下の点に気をつける必要があります • 入る値は必ず整数である • データは１個ではなく105個の数字でやってくる 【ポイント１】 data { int N; データの長さを別の変数として与えると 何かと都合が良いことが多い array[N] int X; } data { 【ポイント２】 変数を配列（array）で与える場合 array[配列の長さ] (変数型) (変数名); の順で書く必要がある 05 基本的なベイズ推論(2) array[105] int X; } 今回の場合これでも良いのだが • データのサイズが変わる • 同じ長さの別の変数が 今後追加される といったケースを想定するため 左のような書き方が一般的 38
stanコードを書いていこう data { どんな形のデータ(𝑌)が与えられるかを指定する。今回の例では int N; • 各年に起きた大地震の数（105年分） array[N] int X; が与えられている。 ※stanコードでは「こんな形式のデータが来る」を指定します。 実際のデータはRから渡します。 } A parameters { 推定するパラメータ(𝜃)を指定する。今回の例では • 平均（ポアソン分布の𝜆） } A model { の1つだけ。 実際に事後分布の形を規定するもの(𝑃 𝑌 𝜃 𝑃 𝜃 )を指定する。 すなわち事前分布𝑃 𝜃 と尤度𝑃 𝑌 𝜃 の両方を書くブロック。 } 𝑃 𝜃𝑌 = 𝑃 𝑌𝜃 𝑃 𝜃 𝑃 𝑌 05 基本的なベイズ推論(2) 39
parametersブロック 今回推定したいのは「回数の平均」（ポアソン分布のパラメータ） 以下の点に気をつける必要があります • 入る値は必ず整数とは限らない • でも必ず0以上の値をとる parameters { real < A lower=0 > lambda; } 今回は下限だけを指定したらよい 05 基本的なベイズ推論(2) 40
stanコードを書いていこう data { どんな形のデータ(𝑌)が与えられるかを指定する。今回の例では int N; • 各年に起きた大地震の数（105年分） array[N] int X; が与えられている。 ※stanコードでは「こんな形式のデータが来る」を指定します。 実際のデータはRから渡します。 } parameters { real <lower=0> lambda; } A model { 推定するパラメータ(𝜃)を指定する。今回の例では • 平均（ポアソン分布の𝜆） の1つだけ。 実際に事後分布の形を規定するもの(𝑃 𝑌 𝜃 𝑃 𝜃 )を指定する。 すなわち事前分布𝑃 𝜃 と尤度𝑃 𝑌 𝜃 の両方を書くブロック。 } 𝑃 𝜃𝑌 = 𝑃 𝑌𝜃 𝑃 𝜃 𝑃 𝑌 05 基本的なベイズ推論(2) 41
modelブロック 今回の尤度𝑃 𝑌 𝜃 はポアソン分布、事前分布𝑃 𝜃 はガンマ分布 尤度はすべてのデータについて計算してかける必要がある 𝐿(𝜆|𝒙) = 𝑃 𝒙 𝜆 = 𝑃 𝑥1 𝜆 × 𝑃 𝑥2 𝜆 × ⋯ × 𝑃(𝑥105 |𝜆) 【事前分布】 ほぼ無情報事前分布 （データ0.01個分の自信） 【尤度】 配列の何番目かを指定 するときも[ ]を使う model { lambda ~ gamma(0.1,0.01); X[1] ~ poisson(lambda); X[2] ~ poisson(lambda); X[3] ~ poisson(lambda); X[4] ~ poisson(lambda); ︙ X[104] ~ poisson(lambda); X[105] ~ poisson(lambda); } 05 基本的なベイズ推論(2) 各年の発生件数が 独立にポアソン分布に従うので… 42
配列に対する簡単な書き方 方法１：forループを使う model { lambda ~ gamma(0.1,0.01); for (i in 1:105) { （forループのための一時変数）iの値を 1:105つまり … と変えながら {}の中の処理を順番に行う X[i] ~ poisson(lambda); } やっていることとしては X[1] ~ poisson(lambda); X[2] ~ poisson(lambda); X[3] ~ poisson(lambda); X[4] ~ poisson(lambda); ︙ X[104] ~ poisson(lambda); X[105] ~ poisson(lambda); ここがi } 05 基本的なベイズ推論(2) 43
配列に対する簡単な書き方 方法２：ベクトル化された関数を使う 実は今回の場合，以下の書き方でOK stanのmodelブロックでは，上から順に model { 対数尤度を計算して足していっている lambda ~ gamma(0.1,0.01); 【例】データが２個だけの場合 X ~ poisson(lambda); } log 𝑃(𝑥1 |𝜆) + log 𝑃(𝑥2 |𝜆) + ⋯ + log 𝑃(𝑥𝑛 |𝜆) わかりやすい上に 計算時間も短くて済む書き方です 𝑃 𝜆𝒙 ∝𝑃 𝒙𝜆 𝑃 𝜆 = 𝑃 𝑥1 𝜆 𝑃 𝑥2 𝜆 ⋯ 𝑃 𝑥𝑛 𝜆 𝑃 𝜆 model { lambda ~ gamma(0.1,0.01); X[1] ~ poisson(lambda); X[2] ~ poisson(lambda); } log 𝑃 𝜆 + log 𝑃(𝑥1 |𝜆) + log 𝑃(𝑥2 |𝜆) stanに実装されている確率分布の関数のほとんどは データが配列で与えられると 対数尤度を自動的にすべて足してくれる 05 基本的なベイズ推論(2) 44
stanコードの完成
data {
どんな形のデータ(𝑌)が与えられるかを指定する。今回の例では
int N;
• 各年に起きた大地震の数（105年分）
array[N] int X;
が与えられている。
※stanコードでは「こんな形式のデータが来る」を指定します。
実際のデータはRから渡します。
}
parameters {
real <lower=0> lambda;
• 平均（ポアソン分布の𝜆）
の1つだけ。
}
model {
lambda ~ gamma(0.1,0.01);
X ~ poisson(lambda);
}
推定するパラメータ(𝜃)を指定する。今回の例では
実際に事後分布の形を規定するもの(𝑃 𝑌 𝜃 𝑃 𝜃 )を指定する。
すなわち事前分布𝑃 𝜃 と尤度𝑃 𝑌 𝜃 の両方を書くブロック。
𝑃 𝜃𝑌 =
𝑃 𝑌𝜃 𝑃 𝜃
𝑃 𝑌
05 基本的なベイズ推論(2)
45
完成したモデル
model_poisson.stan
data {
int N;
array[N] int X;
}
推定に必要な情報
データ
𝑌
各年におきた 𝑥 = (0,1,1,3,21,3, ⋯ )
地震の数
推定したい
パラメータ
𝜃
年間の平均
地震発生件数
尤度
parameters {
real <lower=0> lambda;
今回の事例
𝑃(𝑌|𝜃) ポアソン分布
事前分布
ガンマ分布
𝑃(𝜃)
𝜆
}
𝐺𝑎𝑚𝑚𝑎(𝜆|𝛼, 𝛽)
▼を簡略化した書き方です
𝑥𝑖
𝑖 data
lambda ~ gamma(0.1,0.01);
X ~ poisson(lambda);
𝑃𝑜𝑖𝑠(𝑥|𝜆)
【今回のベイズモデリングのplate notation】
}
model {
𝜆
パラメータ𝜆は
全データで共通
05 基本的なベイズ推論(2)
𝑥1
𝜆
𝑥2
︙
𝑥105
46
あとは実行するだけ
model_poisson.stan
data {
int N;
array[N] int X;
}
推定に必要な情報
データ
𝑌
各年におきた 𝑥 = (0,1,1,3,21,3, ⋯ )
地震の数
推定したい
パラメータ
𝜃
年間の平均
地震発生件数
尤度
parameters {
今回の事例
事前分布
𝜆
𝑃(𝑌|𝜃) ポアソン分布
𝑃(𝜃)
ガンマ分布
𝑃𝑜𝑖𝑠(𝑥|𝜆)
𝐺𝑎𝑚𝑚𝑎(𝜆|𝛼, 𝛽)
real <lower=0> lambda;
}
library(cmdstanr)
model {
model <- cmdstan_model("model_poisson.stan")
lambda ~ gamma(0.1,0.01);
stan_data <- list(N=105, X=dat$time)
X ~ poisson(lambda);
result <- model$sample(data = stan_data)
}
05 基本的なベイズ推論(2)
47
結果を見ていく 【事後分布の要約】 その他の点推定・区間推定の方法は result$summary() 資料04で確認してください 対数尤度 左から 事後平均値(EAP)，事後中央値（MED），事後分布のSD 90%確信区間(Equal-tailed interval) 【事後分布のプロット】 完 全 に 一 致 result$draws() |> mcmc_dens(pars = "lambda") 05 基本的なベイズ推論(2) 𝐺𝑎𝑚𝑚𝑎 590.1, 105.01 48
２ 分析実践編(3) 正規分布のパラメータ 05 基本的なベイズ推論(2) 49
事例 例 あるコンビニチェーンのアナリストは，各店舗の利益に影響する要因を調べることにしました。 その第一歩として，まずは母集団（全店舗）での利益の平均と分散を推測したいと思います。 なお，利益はふつう正規分布に従うと言われているとします。 【データの読み込み】 この名前はなんでもOK “data_cvs.csv”をワーキングディレクトリに配置して 中身はこんなデータ dat <- read.csv("data_cvs.csv") 今日はこれだけ使います sales: その店の一日あたり平均利益（単位：千円） dist: 最寄り駅からの距離（単位：km） floor: 床面積（単位：m2） items: 取扱いアイテム数（単位：個） 全国10の地域から region: その店舗のある地域 それぞれ10店舗ずつ無作為抽出 neighbor: 半径1km以内にあるコンビニの数 したという想定です ※データは適当に作ったので、実際とは異なります。 05 基本的なベイズ推論(2) 50
事例 例 あるコンビニチェーンのアナリストは，各店舗の利益に影響する要因を調べることにしました。 その第一歩として，まずは母集団（全店舗）での利益の平均と分散を推測したいと思います。 なお，利益はふつう正規分布に従うと言われているとします。 【まずは事例の整理】 推定に必要な情報 データ 今回の事例 𝑌 推定したいパラメータ 𝜃 尤度 𝑃(𝑌|𝜃) 事前分布 𝑃(𝜃) 各店舗の利益 𝑥 = (6.52,8.53, ⋯ ) 利益の平均 𝜇 利益の分散 𝜎2 正規分布 𝑁𝑜𝑟𝑚𝑎𝑙(𝑥|𝜇, 𝜎 2 ) 後ほど紹介します 05 基本的なベイズ推論(2) 51
念のため正規性の確認 例 あるコンビニチェーンのアナリストは，各店舗の利益に影響する要因を調べることにしました。 その第一歩として，まずは母集団（全店舗）での利益の平均と分散を推測したいと思います。 なお，利益はふつう正規分布に従うと言われているとします。 【ちゃんとデータを確認する】 そもそも今回のデータがきちんと正規分布に従っていると言えないと 正規分布のパラメータを推定するのはおかしい もしそもそも正規分布に従っていない場合 hist(dat$sales) • サンプリングが偏っていたのか？ • そもそも正規分布に従わないものなのか？ 考えた上で分析モデルを修正する必要がある 今回は正規分布に従っていたとみなします 05 基本的なベイズ推論(2) 52
（補足）正規性の確認 ヒストグラムはガタガタしていると思うあなたへ plot(density(dat$sales)) データがそれぞれ独立して正規分布に従うかどうかの検定 コルモゴロフ・スミノルフ検定 ks.test(dat$sales,y="pnorm",mean=mean(dat$sales),sd=sd(dat$sales)) シャピロ・ウィルク検定 shapiro.test(dat$sales) いずれも帰無仮説を「正規分布に従う」とした検定を行います 05 基本的なベイズ推論(2) 53
非ベイズ的点推定（最尤推定） 正規分布の尤度＝確率密度関数 𝑥−𝜇 2 𝐿 𝜇, 𝜎 𝑥 = 𝑃 𝑥 𝜇, 𝜎 = exp − 2 2𝜎 2 2𝜋𝜎 1 尤度関数は，パラメータが複数あるときには多変数関数として扱います 各店舗の利益が独立であると仮定すると 𝑛 資料03 pp. 32-33 等高線プロット 𝑥𝑖 − 𝜇 2 𝐿 𝜇, 𝜎 𝒙 = ෑ exp − 2 2𝜎 2 2𝜋𝜎 𝑖=1 1 この関数が最大値を取る 𝜇, 𝜎 2 の値を探します 05 基本的なベイズ推論(2) 54
非ベイズ的点推定 𝑛 対数を取って 𝑛 1 𝐿𝐿(𝜇, 𝜎|𝒙) = log 2𝜋𝜎 2 −2 + 𝑖=1 𝑖=1 𝑥𝑖 − 𝜇 2 − 2𝜎 2 𝑛 𝑛 1 2 = − log 2𝜋𝜎 − 2 𝑥𝑖 − 𝜇 2 2 2𝜎 𝑖=1 𝜇 と 𝜎 でそれぞれ偏微分すると 𝑛 • 𝜇 で偏微分：第1項はまるごと消えて，最終的に 𝜕𝐿𝐿(𝜇, 𝜎|𝒙) 1 = 2 𝑥𝑖 − 𝜇 𝜕𝜇 𝜎 𝑖=1 𝑛 𝜕𝐿𝐿(𝜇, 𝜎|𝒙) 𝑁 1 2 = − + 𝑥 − 𝜇 • 𝜎 2 で偏微分： 𝑖 𝜕𝜎 2 2𝜎 2 2 𝜎 2 2 𝑖=1 05 基本的なベイズ推論(2) 55
非ベイズ的点推定 これらが同時に0となるところが（対数）尤度関数の最大なので 【最尤推定値】 𝑛 1 𝑥𝑖 − 𝜇Ƹ = 0 2 𝜎 𝑛 1 𝜇Ƹ = 𝑥𝑖 𝑛 𝑖=1 𝑁 𝑁 1 2 =0 − 2+ 𝑥 − 𝜇 Ƹ 𝑖 2𝜎 2 𝜎2 2 𝑛=1 つまり標本平均 𝑖=1 𝑛 𝜎ො 2 = 1 𝑥𝑖 − 𝜇Ƹ 2 𝑛 つまり標本分散 𝑖=1 05 基本的なベイズ推論(2) 56
偏微分補足 𝑛 𝑛 1 2 𝐿𝐿(𝜇, 𝜎|𝒙) = − log 2𝜋𝜎 − 2 𝑥𝑖 − 𝜇 2 2 2𝜎 𝑖=1 まずは 𝜇 で偏微分 第1項には𝜇がないので無視する 第2項について見てみると 𝑛 1 1 2 − 2 𝑥𝑖 − 𝜇 = − 2 2𝜎 2𝜎 𝑥1 − 𝜇 2 + 𝑥2 − 𝜇 2 + ⋯ + 𝑥𝑛 − 𝜇 2 𝑖=1 𝑓 𝑥 𝑛 を微分すると𝑛𝑓 𝑥 𝑛−1 𝑓′(𝑥)になるので， 𝑥𝑖 − 𝜇 2 を微分すると−2 𝑥𝑖 − 𝜇 1 = − 2 −2 𝑥1 − 𝜇 − 2 𝑥2 − 𝜇 − ⋯ − 2 𝑥𝑛 − 𝜇 2𝜎 05 基本的なベイズ推論(2) 𝑛 1 = 2 𝑦𝑖 − 𝜇 𝜎 𝑖=1 57
偏微分補足 𝑛 𝑛 1 2 𝐿𝐿(𝜇, 𝜎|𝒙) = − log 2𝜋𝜎 − 2 𝑥𝑖 − 𝜇 2 2 2𝜎 𝑖=1 つぎは 𝜎 2 で偏微分 𝑛 𝑁 𝑁 1 2 𝐿𝐿(𝜇, 𝜎|𝒙) = − log 2𝜋 − log 𝜎 − 2 𝑥𝑖 − 𝜇 2 2 2 2𝜎 𝑖=1 消える 𝑛 𝜕𝐿𝐿(𝜇, 𝜎|𝒙) 𝑁 1 2 = − − 𝑥 − 𝜇 𝑖 𝜕𝜎 2 2𝜎 2 2 𝜎 2 2 𝑖=1 𝑑log(𝑥) 1 = 𝑑𝑥 𝑥 1 = 𝑥 −1 なので微分すると 𝑥 1 −1𝑥 −2 = − 2 𝑥 05 基本的なベイズ推論(2) 58
非ベイズ的点推定 ということでデータの平均と分散を求めてみる mean(dat$sales) var(dat$sales) [1] 6.9842 [1] 1.684208 一店舗の利益の確率分布 𝑁 𝑥 𝜇 = 6.98, 𝜎 = 1.68 var()が返す値は不偏分散です。 最尤推定値は標本分散なので補正が必要になります。 【問】 この正規分布に従うとき利益が5千円未満の店舗の割合は 【答】 pnorm(5, 6.98, sqrt(1.68)) [1] 0.06330555 およそ6.33% 05 基本的なベイズ推論(2) 59
非ベイズ的区間推定①平均パラメータ 1 とりあえず95%区間を作る 求めたい区間の上限・下限をそれぞれ𝜇𝐿 , 𝜇𝑈 とする 𝑃 𝜇𝐿 ≤ 𝜇 ≤ 𝜇𝑈 = 0.95 になるような𝜇𝐿 , 𝜇𝑈 の値を求めたら良い 𝜇𝐿 , 𝜇𝑈 の値をどのように設定したら「𝜇𝐿 から𝜇𝑈 の間に母平均 𝜇 が含まれている確率（割合）が95%になる」のかを求めたい 2 既知の確率分布に従う統計量になるように変形する 𝜎 母分散が既知の場合，中心極限定理により，標本平均の標本分布は𝑋ത ∼ 𝑁 𝜇, 2 𝑛 母分散が未知の場合，代わりに不偏推定量である不偏分散𝑠ǁ𝑥2 = となるが， 𝑛 𝑠𝑥2 を使う必要がある 𝑛−1 ത 𝑋−𝜇 標本平均を不偏分散を用いて標準化した𝑍ҧ = は，自由度𝑛 − 1の 𝑡 分布に従う 𝑍ҧ = ǁ 𝑠𝑥 / 𝑛 𝑋ത − 𝜇 ∼ 𝑡(𝑛 − 1) 𝑠ǁ𝑥 / 𝑛 一旦逆になりますが気にしない 𝑃 𝜇𝐿 ≤ 𝜇 ≤ 𝜇𝑈 ത 𝑋−𝜇 の真ん中が ǁ になるように変形させると 𝑠𝑥 / 𝑛 𝑃 05 基本的なベイズ推論(2) ത 𝑈 ത ത 𝐿 𝑋−𝜇 𝑋−𝜇 𝑋−𝜇 ≤ ǁ ≤ ǁ 𝑠ǁ𝑥 / 𝑛 𝑠𝑥 / 𝑛 𝑠𝑥 / 𝑛 60
非ベイズ的区間推定①平均パラメータ 3 もう一つ95%区間を作る ത 𝑋−𝜇 𝑍ҧ = ǁ が自由度𝑛 − 1の 𝑡 分布に従うということは（今回のサンプルサイズは100なので） 𝑠𝑥 / 𝑛 ത 𝑋−𝜇 ≤ 1.984 𝑠𝑥 / 𝑛 𝑡分布表を用いると，𝑃 −1.984 ≤ ǁ = 0.95 と分かる 4 2つの式を対応させると… 𝑋ത − 𝜇 ≤ 1.96 = 0.95 3 より 𝑃 −1.96 ≤ 𝑠ǁ𝑥 / 𝑛 𝑋ത − 𝜇𝑈 𝑋ത − 𝜇 𝑋ത − 𝜇𝐿 ≤ ≤ = 0.95 2 より 𝑃 𝑠ǁ𝑥 / 𝑛 𝑠ǁ𝑥 / 𝑛 𝑠ǁ𝑥 / 𝑛 𝑠ǁ𝑥 𝑠ǁ𝑥 𝑃 𝑋ത − 1.984 ≤ 𝜇 ≤ 𝑋ത + 1.984 = 0.95 𝑛 𝑛 すべての標本でこの区間を作った場合 95%の割合で真値𝜇が含まれる 𝑋ത − 𝜇𝐿 = 1.984 𝑠ǁ𝑥 / 𝑛 𝑋ത − 𝜇𝑈 = −1.984 𝑠ǁ𝑥 / 𝑛 𝑥ҧ = 6.98 𝑠ǁ𝑥 = 1.684 𝑛 = 100 𝜇𝐿 = 𝑋ത − 1.984 𝑠ǁ𝑥 𝑛 𝜇𝑈 = 𝑋ത + 1.984 𝑠ǁ𝑥 𝑛 6.98 − 1.984 × 0.13 ≤ 𝜇 ≤ 6.98 + 1.984 × 0.13 既知の値を 当てはめると 05 基本的なベイズ推論(2) 【答】 およそ6.722から7.238 61
非ベイズ的区間推定②分散パラメータ する 1 とりあえず95%区間を作る = 求めたい区間の上限・下限をそれぞれ𝜎𝐿2 , 𝜎𝑈2 とする 𝑃 𝜎𝐿2 ≤ 𝜎 2 ≤ 𝜎𝑈2 = 0.95 になるような𝜎𝐿2 , 𝜎𝑈2 の値を求めたら良い 平均は 定される は𝑈2 の値をどのように設定したら に関して標準化されたもの 𝜎𝐿2 , 𝜎 「𝜎𝐿2 から 𝜎𝑈2 の間に母分散 𝜎 2 が含まれている確率（割合）が95%になる」のかを求めたい が変われば も変わるため 2 既知の確率分布に従う統計量になるように変形する 1 ) でないといけない 詳しく はこ ちら 𝑥𝑖 を標準化した値を 𝑧𝑖 とすると，𝑧𝑖2 の和は自由度𝑛 − 1 の𝜒 2 分布に従う 2 の和 の分散 2 1 ( 2 1) ( 1) 詳しくはこちら 𝑛 倍したもの 𝑛−1 𝑥 の不偏分散 𝑠ǁ𝑥2 は，標本分散を 𝜎2 2 2 𝑠ǁ𝑥 は 𝜒 (𝑛 − 1) に従う 𝑛−1 𝑃 の分散 標本分散 2 2 ( 1) 𝑛−1 2 2 (𝑛 − 1) に従う 𝑠 ǁ は 𝜒 𝑥 𝜎2 05 基本的なベイズ推論(2) 𝑃 𝜎𝐿2 ≤ 𝜎 2 ≤ 𝜎𝑈2 𝑛−1 2 𝑛−1 2 𝑛−1 2 𝑠ǁ𝑥 ≤ 𝑠ǁ ≤ 𝑠ǁ𝑥 𝜎2 𝑥 𝜎𝑈2 𝜎𝐿2 62
非ベイズ的区間推定②分散パラメータ 3 もう一つ95%区間を作る 𝑛−1 2 𝑠Ƹ が95%の確率で含まれる区間は 𝜎2 𝑥 𝑃 73.36 ≤ 確率密度 𝜒 2 (𝑛 − 1) に従う 𝑛−1 2 𝑠ǁ ≤ 128.42 = 0.95 𝜎2 𝑥 𝜒 2 (99) 今回のサンプルサイズは100なので 自由度99の𝜒 2 分布の上下2.5%点の位置を求める 4 2つの式を対応させると… 3 より 𝑛−1 2 𝑃 73.36 ≤ 𝑠ǁ ≤ 128.42 = 0.95 𝜎2 𝑥 𝑛−1 2 𝑛−1 2 𝑛−1 2 より 𝑃 𝑠ǁ𝑥 ≤ 𝑠ǁ ≤ 𝑠ǁ𝑥 = 0.95 2 𝜎2 𝑥 𝜎𝑈2 𝜎𝐿2 𝑛−1 2 𝑛−1 2 𝑃 𝑠ǁ𝑥 ≤ 𝜎 2 ≤ 𝑠ǁ = 0.95 128.42 73.36 𝑥 すべての標本でこの区間を作った場合 95%の割合で真値𝜎 2 が含まれる 𝑠ǁ𝑥2 = 1.684 𝑛 = 100 𝑛−1 2 𝑠ǁ𝑥 = 128.42 𝜎𝐿2 𝑛−1 2 𝑠ǁ𝑥 = 73.36 𝜎𝑈2 既知の値を 当てはめると 𝜎𝐿2 = 𝑛−1 2 𝑠ǁ 128.42 𝑥 𝜎𝑈2 = 𝑛−1 2 𝑠ǁ 73.36 𝑥 99 99 × 1.684 ≤ 𝜎 2 ≤ × 1.684 128.42 73.36 【答】 およそ1.298から2.273 05 基本的なベイズ推論(2) 63
ベイズ推定（まずは解析的に） パラメータが複数あるけど… 推定に必要な情報 データ 今回の事例 𝑌 推定したいパラメータ 𝜃 尤度 𝑃(𝑌|𝜃) 事前分布 𝑃(𝜃) 各店舗の利益 𝑥 = (6.52,8.53, ⋯ ) 利益の平均 𝜇 利益の分散 𝜎2 正規分布 𝑁𝑜𝑟𝑚𝑎𝑙(𝑥|𝜇, 𝜎 2 ) この場合の事前分布は？ 各店舗の利益が独立であると仮定すると，データ全体の尤度は 𝑛 𝑥𝑖 − 𝜇 2 𝐿 𝜇, 𝜎 𝒙 = ෑ exp − 2 2𝜎 2 2𝜋𝜎 𝑖=1 1 05 基本的なベイズ推論(2) 64
条件付き確率で考えていこう パラメータが複数ある場合，条件付き確率の積に分解して考えることが多い 𝑃 𝒙 𝜇, 𝜎 𝑃(𝜇, 𝜎) 𝑃 𝜇, 𝜎 𝒙 = ∝ 𝑃 𝒙 𝜇, 𝜎 𝑃 𝜇, 𝜎 𝑃(𝒙) = 𝑃 𝒙 𝜇, 𝜎 𝑃 𝜇|𝜎 𝑃 𝜎 この場合，共役事前分布には • 𝜎 が既知の場合の 𝜇 の共役事前分布 • 𝜎 の共役事前分布 の積を用意してあげたら良い！ ということで，まずは 𝜎 が既知の場合の 𝜇 の共役事前分布 を見ていきます 05 基本的なベイズ推論(2) 65
母分散が既知の場合の𝜇の共役事前分布？ 尤度は正規分布｜𝑥𝑖 ∼ 𝑁(𝜇, 𝜎) 𝑛 𝑥𝑖 − 𝜇 2 𝐿 𝜇 𝒙, 𝜎 = 𝑃 𝒙 𝜇, 𝜎 = ෑ exp − 2 2𝜎 2 2𝜋𝜎 𝑖=1 1 この場合，共役事前分布は正規分布になります｜𝜇 ∼ 𝑁(𝜇0 , 𝜎0 ) 𝑃 𝜇 𝜇0 , 𝜎0 = 1 𝜇 − 𝜇0 exp − 2 2𝜎 0 2 2 データを与える前＝スタート時点での 分布のパラメータということで 添字0をつけて表現することが多い 2𝜋𝜎0 05 基本的なベイズ推論(2) 66
（補足）事後分布の導出 ここから数枚の補足スライドによって 母平均パラメータ𝜇の事後分布が正規分布𝑁 𝜇∗ , 𝜎 ∗ になること を示していきます｜𝑃(𝜇|𝒙) ∼ 𝑁 𝜇∗ , 𝜎 ∗ 事後分布∝事前分布×尤度 事後分布 事前分布 𝑃 𝜇 𝒙, 𝜇0 , 𝜎0 ∝ 1 𝑛 尤度 𝜇 − 𝜇0 2 1 𝑥𝑖 − 𝜇 2 exp − ×ෑ exp − 2 2 2 2𝜎 2𝜎 2𝜋𝜎 0 2 𝑖=1 2𝜋𝜎0 正規化定数 = 𝑛 1 1 2𝜋𝜎02 2𝜋𝜎 2 𝜇 − 𝜇0 2 σ𝑛𝑖=1 𝑥𝑖 − 𝜇 2 × exp − − 2 2𝜎 2 2𝜎0 𝜇 − 𝜇0 2 σ𝑛𝑖=1 𝑥𝑖 − 𝜇 2 ∝ exp − − 2 2𝜎 2 2𝜎0 05 基本的なベイズ推論(2) 67
（補足）事後分布の導出 𝑃 𝜇 𝒙, 𝜇0 , 𝜎0 𝜇 − 𝜇0 2 σ𝑛𝑖=1 𝑥𝑖 − 𝜇 2 ∝ exp − − 2 2𝜎 2 2𝜎0 指数の中身を𝜇について整理していく 𝜇 − 𝜇0 2 σ𝑛𝑖=1 𝑥𝑖 − 𝜇 2 1 − − = − 2𝜎 2 2 2𝜎02 𝑃 𝜇 𝒙, 𝜇0 , 𝜎0 𝜇 と無関係 𝑛 𝑛 2 2 σ σ 1 𝑛 𝜇 𝑥 𝜇 𝑥 0 𝑖 0 𝑖 𝑖=1 𝑖=1 2 𝜇+ 2+ 2 + 𝜎2 𝜇 + 2 2 + 2 𝜎2 𝜎2 𝜎0 𝜎0 𝜎0 1 ∝ exp − 2 𝑛 σ 1 𝑛 𝜇 0 𝑖=1 𝑥𝑖 2 𝜇 2 + 𝜎2 𝜇 + 2 2+ 2 𝜎 𝜎0 𝜎0 最終的に正規分布のカーネルと同じ形 1 𝜇 − 𝜇∗ 2 exp − 2 𝜎∗2 になるように，さらに変形していきます 05 基本的なベイズ推論(2) 68
（補足）事後分布の導出 𝑃 𝜇 𝒙, 𝜇0 , 𝜎0 ここで，𝑎 = 1 𝑛 + 𝜎02 𝜎2 1 ∝ exp − 2 ，𝑏 = 𝑛 σ 1 𝑛 𝜇 0 𝑖=1 𝑥𝑖 2 + 𝜇 + 2 + 𝜇 𝜎2 𝜎02 𝜎 2 𝜎02 σ𝑛 𝜇0 𝑖=1 𝑥𝑖 + 𝜎02 𝜎2 𝑃 𝜇 𝒙, 𝜇0 , 𝜎0 ∝ exp − 【最終目標】 1 𝜇 − 𝜇∗ 2 𝑃 𝜇 𝒙, 𝜇0 , 𝜎0 ∝ exp − 2 𝜎∗2 とおく 1 𝑎𝜇2 + 2𝑏𝜇 2 1 𝑏 = exp − 𝑎 𝜇 − 2 𝑎 1 𝑏 ∝ exp − 𝑎 𝜇 − 2 𝑎 2 2 𝜇 と無関係 1 𝑏2 + 2𝑎 𝑏 2 1 𝜇−𝑎 = exp − 2 1 2 𝑎 05 基本的なベイズ推論(2) 69
母分散が既知の場合の𝜇の共役事前分布 共役事前分布として正規分布｜𝜇 ∼ 𝑁(𝜇0 , 𝜎0 )をおいた場合の事後分布は 𝑃 𝜇 𝒙, 𝜇0 , 𝜎0 𝑎= 1 𝑛 + 2 𝜎02 𝜎 𝑏 2 1 𝜇−𝑎 ∝ exp − 2 1 2 𝑎 ，𝑏 = σ𝑛 𝜇0 𝑖=1 𝑥𝑖 + 2 𝜎02 𝜎 𝜇0 𝑛𝑥ҧ 2 + 𝜎2 𝜎0 𝜇𝑝𝑜𝑠𝑡 ∼ 𝑁 , 1 𝑛 + 𝜎02 𝜎 2 𝑏 𝜇𝑝𝑜𝑠𝑡 ∼ 𝑁 , 𝑎 𝑎 としていたので 𝑛 𝑥𝑖 = 𝑛𝑥ҧ 1 𝑖=1 1 𝑛 + 𝜎02 𝜎 2 05 基本的なベイズ推論(2) 見方を変えてみましょう 70
分散の逆数 （統計的・ベイズ）推定の文脈における確率分布の分散 予測（や信念）の精度を表している，と見ることができる 予測の精度があまり高くない 𝑁(0, 2) -3かもしれないし 4とかかもしれない 予測の精度が結構高い 𝑁(0, 0.5) ほぼ確実に -1から1の間 分散の逆数を「精度」と呼ぶことがあります 05 基本的なベイズ推論(2) 71
母分散が既知の場合の𝜇の共役事前分布 𝜇0 𝑛𝑥ҧ 2 + 𝜎2 𝜎0 𝜇𝑝𝑜𝑠𝑡 ∼ 𝑁 , 1 𝑛 2 + 𝜎2 𝜎0 𝑃 𝜇 = 𝑁(𝜇0 , 𝜎0 ) 1 1 𝑛 + 𝜎02 𝜎 2 「事前に持っていた信念の精度」を𝜏0 = 1 「母集団分布の精度」を𝜏 = 2 𝜎 𝜏0 𝜇0 + 𝑛𝜏𝑥ҧ 1 𝜇𝑝𝑜𝑠𝑡 ∼ 𝑁 , 𝜏0 + 𝑛𝜏 𝜏0 + 𝑛𝜏 1 𝜎02 …とおくと 見覚えのある形？ つまり 𝑁 𝜇0 , 1 𝜏0 が尤度によって更新され𝑁 05 基本的なベイズ推論(2) 𝜏0 𝜇0 +𝑛𝜏𝑥ҧ 1 , になった 𝜏0 +𝑛𝜏 𝜏0 +𝑛𝜏 72
正規分布の更新とパラメータの解釈 𝑁 𝜇0 , 例 1 𝜏0 𝑁 10, 𝜏0 𝜇0 +𝑛𝜏𝑥ҧ 1 , になった 𝜏0 +𝑛𝜏 𝜏0 +𝑛𝜏 が尤度によって更新され𝑁 1 𝑁 0.01 0.01 × 10 + 100𝜏 × 6.98 1 , 0.01 + 100𝜏 0.01 + 100𝜏 × 6.98 • かなり精度𝜏0 が低い信念 事後平均は事前とデータの重み付け和 • 平均は10くらいだろう 事後分散は事前とデータの精度の和 データ 尤度 100店舗で 𝑥ҧ = 6.98 05 基本的なベイズ推論(2) 73
信念の強さの比較 【事前の信念】 予想 【データ】 100店舗の平均が6.98 平均10くらいだろう 正規分布の 𝜇 = 6.98 と予想 正規分布の 𝜇 = 10 と予想 自信 順当に行けば 100店舗分のデータがあるので 「正直言って全く自信はないです」 精度0.01（ほぼゼロ） 精度100𝜏 【例】 母分散が2の場合 𝜏 = 0.5になるので データが与える精度は50 事前の信念はほぼ自信なし（無情報に近い事前分布）だったので 推論の結果はほぼ完全にデータ（尤度）によって決定されました。 05 基本的なベイズ推論(2) 74
更新前後の期待値 事前分布 𝑁 𝜇0 , 1 𝜏0 事後分布 が尤度によって更新され𝑁 𝜏0 𝜇0 +𝑛𝜏𝑥ҧ 1 , になった 𝜏0 +𝑛𝜏 𝜏0 +𝑛𝜏 正規分布𝑁 𝜇, 𝜎 の期待値は 𝜇 𝜏0 𝜇0 + 𝑛𝜏𝑥ҧ 𝜏0 𝑛𝜏 = 𝜇0 + 𝑥ҧ 事後分布の期待値は 𝜏0 + 𝑛𝜏 𝜏0 + 𝑛𝜏 𝜏0 + 𝑛𝜏 事前 分布 尤度 パラメータ 意味 𝜏0 事前情報の総量 𝜇0 事前期待値 𝑛𝜏 データの 情報の総量 𝑥ҧ データの平均値 合計1 事後期待値 = 事前情報のウェイト × 事前期待値 + データのウェイト × (データの平均値) 05 基本的なベイズ推論(2) 75
ようやく1個終わりました パラメータが複数ある場合，条件付き確率の積に分解して考えることが多い 𝑃 𝒙 𝜇, 𝜎 𝑃(𝜇, 𝜎) 𝑃 𝜇, 𝜎 𝒙 = ∝ 𝑃 𝒙 𝜇, 𝜎 𝑃 𝜇, 𝜎 𝑃(𝒙) = 𝑃 𝒙 𝜇, 𝜎 𝑃 𝜇|𝜎 𝑃 𝜎 この場合，共役事前分布には 正規分布 1 𝑁 𝜇0 , = 𝑁(𝜇0 , 𝜎0 ) 𝜏0 • 𝜎 が既知の場合の 𝜇 の共役事前分布 • 𝜎 の共役事前分布 の積を用意してあげたら良い！ 続いて 𝜎 の共役事前分布 を見ていきます 05 基本的なベイズ推論(2) 76
とその前に まずは𝜇の事前分布を少しだけ書き換えておきます 後ほどパラメータの解釈のときに使うため 𝑃 𝜇 = 𝑁(𝜇0 , 𝜎0 ) 𝜎0 は「事前に持っている信念の強さ（精度）」の逆数 もし事前に大量の情報を持っていた場合， 𝜎0 は小さくなるはず 「事前の情報量（サンプルサイズ）」を表すパラメータ 𝑛0 を使う 標本理論で言うところの サンプルサイズが大きいほど 標本分布の分散が小さくなる仕組み と同じ 母分散 𝜎 については？ 今回は母分散は未知 標本理論で言うところの 中心極限定理において ただし母分散が大きいほど母平均の推定のばらつきは大きくなるはず 分散の分子が母分散であること と同じ 以上を組み合わせて 𝑃 𝜇|𝜎 = 𝑁 𝜇0 , 05 基本的なベイズ推論(2) 𝜎 𝑛0 と表しておきます 77
結局共役事前分布は？ 𝑃 𝜇, 𝜎 = 𝑃 𝜇 𝜎 𝑃(𝜎 2 ) 2つのパラメータ𝜇, 𝜎の 同時事前分布 ＝ 𝜎を所与とした𝜇の事前分布 ＝正規分布 𝑃 𝜇|𝜎 = 𝑁 𝜇0 , 𝜎 𝑛0 パラメータ𝜎 2 の事前分布 × ＝逆ガンマ分布 𝑃 𝜎 2 𝜈0 𝜈0 𝜎0 = 𝐼𝑛𝑣𝐺𝑎𝑚𝑚𝑎 , 2 2 突如現れた「逆ガンマ分布」とは何者か？ そしてパラメータ𝜈0 , 𝜎0 の意味は？ 05 基本的なベイズ推論(2) 78
逆ガンマ分布｜Inverse-Gamma distribution ガンマ分布に従う確率変数の逆数が従う分布 ガンマ分布の確率密度関数は 𝛽𝛼 𝑃 𝑥 𝛼, 𝛽 = Γ 𝛼 関数 パラメータ 𝛼 𝛽 1 𝑥 𝛼+1 （よくわからない） 𝛽 𝛼 𝛼−1 −𝛽𝑥 𝑃 𝑥 𝛼, 𝛽 = 𝑥 𝑒 Γ 𝛼 1 −𝛽 𝑒 𝑥 【ここでのポイント】 略記 𝐼𝑛𝑣𝐺𝑎𝑚𝑚𝑎(𝑥|𝛼, 𝛽) 「分散 𝜎 2 が逆ガンマ分布に従う」 期待値 𝛽 𝛼−1 𝛽2 𝛼−2 𝛼−1 2 ということは 分散 つまり 𝑃 𝜎2 「精度（＝分散の逆数）𝜏はガンマ分布に従う」 と言うこともできる 𝜈0 𝜈0 2 𝜈0 𝜈0 𝜎02 = 𝐼𝑛𝑣𝐺𝑎𝑚𝑚𝑎 , 𝜎0 ⟺ 𝑃 𝜏 = 𝐺𝑎𝑚𝑚𝑎 , 2 2 2 2 05 基本的なベイズ推論(2) 79
事後分布の導出①平均パラメータ 先程の導出をほぼそのまま使えば良い 事前分布 p.72 1 𝜏0 𝑁 𝜇0 , 事後分布 が尤度によって更新され𝑁 𝜎 今回は 𝑃 𝜇|𝜎 = 𝑁 𝜇0 , 𝑛0 と置き直しているので… 事前分布 𝑁 𝜇0 , 𝜎 𝑛0 意味 事前 分布 𝑛0 事前情報の総量 𝜇0 事前期待値 尤度 𝑛 データの総量 𝑥ҧ データの平均値 1 𝜏 = 2 であることを利用すると… 𝜎 事後分布 が尤度によって更新され𝑁 パラメータ 𝜏0 𝜇0 +𝑛𝜏𝑥ҧ 1 , になった 𝜏0 +𝑛𝜏 𝜏0 +𝑛𝜏 𝑛0 𝜇0 +𝑛𝑥ҧ 𝜎 , になった 𝑛0 +𝑛 𝑛0 +𝑛 事後期待値 = 事前情報のウェイト × 事前期待値 + データのウェイト × (データの平均値) 05 基本的なベイズ推論(2) 80
事後分布の導出②分散パラメータ 最終的に事後分布は 事前分布 𝐼𝑛𝑣𝐺𝑎𝑚𝑚𝑎 事後分布 𝜈0 𝜈0 𝜎02 , 2 2 ここで が尤度によって更新され𝐼𝑛𝑣𝐺𝑎𝑚𝑚𝑎 𝜈𝑛 𝜈𝑛 𝜎𝑛2 , 2 になる 2 𝜈𝑛 = 𝜈0 + 𝑛 1 𝑛0 𝑛 2 2 2 𝜎𝑛 = 𝜈0 𝜎0 + 𝑛 − 1 𝑠𝑥 + 𝑥ҧ − 𝜇0 2 𝜈𝑛 𝑛0 + 𝑛 よくわからないので精度に変換して考えてみます 事前分布 𝐺𝑎𝑚𝑚𝑎 𝜈0 𝜈0 𝜎02 , 2 2 事後分布 が尤度によって更新され𝐺𝑎𝑚𝑚𝑎 05 基本的なベイズ推論(2) 𝜈𝑛 𝜈𝑛 𝜎𝑛2 , 2 になる 2 81
パラメータの解釈 事前分布 𝐺𝑎𝑚𝑚𝑎 事後分布 𝜈0 𝜈0 𝜎02 , 2 2 が尤度によって更新され𝐺𝑎𝑚𝑚𝑎 ガンマ分布𝐺𝑎𝑚𝑚𝑎(𝛼, 𝛽)のパラメータのポイント 𝛼 分散は 2 である 𝛽 𝜈𝑛 𝜈𝑛 𝜎𝑛2 , 2 になる 2 𝐺𝑎𝑚𝑚𝑎 𝜈𝑛 𝜈𝑛 𝜎𝑛2 , 2 2 p. 12 の分散は𝜈𝑛 に対応して小さくなっていく これを踏まえると… 𝜈𝑛 = 𝜈0 + 𝑛 事後の情報の量 ＝ 事前の情報の量 ＋ データの量 情報の量が多くなるほど推論の精度が高まる ＝事後分布の分散が小さくなる 05 基本的なベイズ推論(2) 82
パラメータの解釈 事前分布 𝐺𝑎𝑚𝑚𝑎 事後分布 𝜈0 𝜈0 𝜎02 , 2 2 が尤度によって更新され𝐺𝑎𝑚𝑚𝑎 ガンマ分布𝐺𝑎𝑚𝑚𝑎(𝛼, 𝛽)のパラメータのポイント 𝛼 期待値は である 𝛽 𝐺𝑎𝑚𝑚𝑎 𝜈𝑛 𝜈𝑛 𝜎𝑛2 , 2 2 𝜈𝑛 𝜈𝑛 𝜎𝑛2 , 2 になる 2 p. 12 の期待値は𝜎2 = 𝜏𝑛 になる 1 𝑛 これを踏まえると… 1 𝑛0 𝑛 𝜈0 𝜎02 + 𝑛 − 1 𝑠𝑥2 + 𝑥ҧ − 𝜇0 2 𝜈𝑛 𝑛0 + 𝑛 𝑛 𝑛 𝜈0 𝜎02 + 𝑛 − 1 𝑠𝑥2 + 0 𝑥ҧ − 𝜇0 2 𝑛0 + 𝑛 = 𝜈0 + 𝑛 − 1 + 1 𝜎𝑛2 = 分散の事後予測 ∝ 分散の事前予測 事前分布とデータで予想していた値が 大きく異なるということは その重み付け和に付随する 精度 𝜏 の値は小さいであろう 平均値の ＋ データによる分散の予測 ＋ 05 基本的なベイズ推論(2) 事前予測とデータのズレ 83
（補足）分散パラメータの事後分布の導出 いつかきちんと書きたい（今は余裕がない） 05 基本的なベイズ推論(2) 84
まとめると 正規-逆ガンマ分布が共役事前分布 𝑃 𝜇, 𝜎 = 𝑃 𝜇 𝜎 𝑃 𝜎 2 𝜎 𝜈0 𝜈0 𝜎02 = 𝑁 𝜇0 , 𝐼𝑛𝑣𝐺𝑎𝑚𝑚𝑎 , 𝑛0 2 2 すると事後分布も正規-逆ガンマ分布になる 𝑃 𝜇, 𝜎|𝒙 ∝ 𝑃 𝜇 𝜎, 𝒙 𝑃 𝜎 2 𝒙 𝑛0 𝜇0 + 𝑛𝑥ҧ 𝜎 𝜈𝑛 𝜈𝑛 𝜎𝑛2 =𝑁 , 𝐼𝑛𝑣𝐺𝑎𝑚𝑚𝑎 , 𝑛0 + 𝑛 2 2 𝑛0 + 𝑛 𝜈𝑛 = 𝜈0 + 𝑛 1 𝑛0 𝑛 2 2 2 𝜎𝑛 = 𝜈0 𝜎0 + 𝑛 − 1 𝑠𝑥 + 𝑥ҧ − 𝜇0 2 𝜈𝑛 𝑛0 + 𝑛 05 基本的なベイズ推論(2) 85
事前分布のパラメータ設定 𝑃 𝜇, 𝜎|𝒙 ∝ 𝑃 𝜇 𝜎, 𝒙 𝑃 𝜎 𝒙 𝑛0 𝜇0 + 𝑛𝑥ҧ 𝜎 𝜈𝑛 𝜈𝑛 𝜎𝑛2 =𝑁 , 𝐼𝑛𝑣𝐺𝑎𝑚𝑚𝑎 , 𝑛0 + 𝑛 2 2 𝑛0 + 𝑛 平均・分散ともに「事前の情報量」を意味するパラメータがある とりあえずそれらを小さめにしておけば事前分布の影響は小さくなっていく 𝜎 𝜈0 𝜈0 𝜎02 𝑃 𝜇, 𝜎 = 𝑃 𝜇 𝜎 = 𝑁 𝜇0 , 𝐼𝑛𝑣𝐺𝑎𝑚𝑚𝑎 , 𝑛0 2 2 標準偏差を大きくする 05 基本的なベイズ推論(2) 両方のパラメータを 小さくする 86
事後分布の導出例 事前分布をかなり無情報にしてみる(𝑛0 = 𝜈0 = 0.001, 𝜇0 = 0, 𝜎0 = 1) データからそれぞれ統計量を計算する 標本平均 𝑥ҧ mean(dat$sales) [1] 6.9842 標本分散𝑠𝑥2 var(dat$sales) [1] 1.684208 本当は不偏分散 𝑛0 𝜇0 + 𝑛𝑥ҧ 𝜎 𝜎 𝑃 𝜇|𝜎, 𝒙 ∝ 𝑁 , ≃ 𝑁 6.9842, 𝑛0 + 𝑛 10 𝑛0 + 𝑛 𝑛0 𝑛 2 2 2 𝜈 𝜎 + 𝑛 − 1 𝑠 + 𝑥 ҧ − 𝜇 0 𝑥 0 0 𝜈0 + 𝑛 𝑛0 + 𝑛 𝑃 𝜎 2 |𝒙 ∝ 𝐼𝑛𝑣𝐺𝑎𝑚𝑚𝑎 , ≃ 𝐼𝑛𝑣𝐺𝑎𝑚𝑚𝑎(50.0005, 83.383) 2 2 05 基本的なベイズ推論(2) 87
分散パラメータの事後分布 事前分布 𝐼𝑛𝑣𝐺𝑎𝑚𝑚𝑎 例 𝜈0 𝜈0 𝜎02 , 2 2 事後分布 が尤度によって更新され𝐼𝑛𝑣𝐺𝑎𝑚𝑚𝑎 0.001 0.001 × 1 𝐼𝑛𝑣𝐺𝑎𝑚𝑚𝑎 , 2 2 𝜈𝑛 𝜈𝑛 𝜎𝑛2 , 2 になる 2 𝐼𝑛𝑣𝐺𝑎𝑚𝑚𝑎(50.0005, 83.383) • データ0.001個分くらいの強さの信念 データの総量は事前とデータの合計 • 分散は1くらいだろう 分散は事前とデータの重み付け和＋α データ 尤度 100店舗で 𝑠𝑥2 = 1.684 05 基本的なベイズ推論(2) 88
事前の信念が強かったら 事前分布 𝐼𝑛𝑣𝐺𝑎𝑚𝑚𝑎 例 𝜈0 𝜈0 𝜎02 , 2 2 事後分布 が尤度によって更新され𝐼𝑛𝑣𝐺𝑎𝑚𝑚𝑎 100 100 × 1 𝐼𝑛𝑣𝐺𝑎𝑚𝑚𝑎 , 2 2 𝜈𝑛 𝜈𝑛 𝜎𝑛2 , 2 になる 2 𝐼𝑛𝑣𝐺𝑎𝑚𝑚𝑎(100, 133.383) • データ100個分くらいの強さの信念 データの総量は事前とデータの合計 • 分散は1くらいだろう 分散は事前とデータの重み付け和＋α データ 尤度 100店舗で 𝑠𝑥2 = 1.684 05 基本的なベイズ推論(2) 89
事前分布による違い ※分散の予測はとりあえず10（精度の予測が0.1）であるとして… ほぼ 無情報事前分布 𝑃 𝜎 2 = 𝐼𝑛𝑣𝐺𝑎𝑚𝑚𝑎 思想の強い事前分布 0.001 0.001 × 1 , 2 2 𝑃 𝜎 2 = 𝐼𝑛𝑣𝐺𝑎𝑚𝑚𝑎 のときの事後分布 100 100 × 1 , 2 2 のときの事後分布 2 𝜎𝑝𝑜𝑠𝑡 ~𝐼𝐺(100, 133.383) 𝑃 𝜎 2 𝑥 = 𝐼𝑛𝑣𝐺𝑎𝑚𝑚𝑎(50.0005, 83.383) データ 尤度 100店舗で 𝑠𝑥2 = 1.684 05 基本的なベイズ推論(2) 90
平均パラメータの事後分布はどうなる？ 平均パラメータ𝜇の条件付き事後分布は 𝑃 𝜇|𝜎, 𝒙 ∝ 𝑁 𝜎 6.9842, 10 𝜎2は 1に近いかもしれないし 2より大きいかもしれない 𝜎 2 = 1だったら 𝑃 𝜇|𝜎, 𝒙 ∝ 𝑁 6.9842, 1 10 𝜎 2 = 2だったら 𝑃 𝜇|𝜎, 𝒙 ∝ 𝑁 6.9842, 2 10 ふつう関心があるのは平均パラメータ𝜇そのものの事後分布 𝑃 𝜇|𝒙 本当は周辺化して求める必要がある ∞ 𝑃 𝜇|𝒙 = න 𝑃 𝜇|𝜎 2 , 𝒙 𝑃 𝜎 2 |𝒙 𝑑𝜎 2 𝜎2 =0 05 基本的なベイズ推論(2) …が，明らかにしんどいので 91
事後分布からの乱数生成 stanとも同じような考え方です 𝑃 𝜇, 𝜎 2 |𝒙 ∝ 𝑃 𝜇 𝜎 2 , 𝒙 𝑃 𝜎 2 𝒙 条件付き分布になっている場合は順番に ① 𝑃 𝜎 2 𝒙 から𝜎 2 の乱数を作る ここから乱数を大量に作れば とうぜん𝜎 2 = 1.6付近の値が多くなり 大体の値は1から2.5付近になる ② いま作った 𝜎 2 を一つずつ用いて， 𝑃 𝜇 𝜎, 𝒙 から乱数を作る 𝜎 2 の乱数を一つずつ 𝑁 6.9842, 𝜎 10 に入れて，正規分布から乱数を一つずつ作れば ∞ 𝑃 𝜇|𝒙 = න 𝑃 𝜇|𝜎 2 , 𝒙 𝑃 𝜎 2 |𝒙 𝑑𝜎 2 𝜎2 =0 05 基本的なベイズ推論(2) 手順①での 𝜎 2 の各値の出現率が 確率的な重み付け𝑃 𝜎 2 |𝒙 を 再現してくれるため 92
事後分布からの乱数生成 𝑃 𝜇, 𝜎 2 |𝒙 ∝ 𝑃 𝜇 𝜎 2 , 𝒙 𝑃 𝜎 2 𝒙 条件付き分布になっている場合は順番に ① 𝑃 𝜎 2 𝒙 から𝜎 2 の乱数を作る n_draw <- 100000 # いくつ乱数を作るか post_sigma <- 1/rgamma(n_draw, 50.0005, 83.383) ガンマ分布に従う乱数（精度）の逆数 ▲ 事後分布 𝑃(𝜎 2 |𝒙) からのサンプリング ② いま作った 𝜎 2 を一つずつ用いて， 𝑃 𝜇 𝜎, 𝒙 から乱数を作る post_mu <- rnorm(n_draw, 6.9842, sqrt(post_sigma)/10) 𝑁 6.9842, 𝜎 に従う乱数 10 ▲ 周辺事後分布 𝑃(𝜇|𝒙) からのサンプリングと同じこと 05 基本的なベイズ推論(2) 93
（補足）Rの乱数生成関数の挙動 rnorm(n, mean, sd)について 分布 できあがり 𝑁(0,1) -0.708 【mean,sdを1個だけ与えた時】 例｜rnorm(n=5, mean=0, sd=1) 𝑁(0,1) 0.164 𝑁(0,1) -1.054 平均mean,標準偏差sdの正規分布から 𝑁(0,1) 0.252 n個乱数を発生 𝑁(0,1) 0.465 【mean,sdをベクトルで与えた時】 例｜rnorm(n=5, mean=0, sd=1:5) 分布 できあがり 𝑁(0,1) 0.535 それぞれ異なるパラメータの正規分布から 𝑁(0,2) -0.577 1個ずつ乱数を発生 𝑁(0,3) 4.320 𝑁(0,4) 2.046 𝑁(0,5) 10.620 05 基本的なベイズ推論(2) SDが大きい 分布から 発生するため 大きな値 94
事後分布のプロット どちらを使ってもOKです ヒストグラム カーネル密度推定 hist(post_sigma) plot(density(post_mu)) 05 基本的なベイズ推論(2) 95
事後分布からの点推定・区間推定 pp.6804 資参 も 料照69 資料04 pp.68- 69も参照 tidybayesパッケージには最頻値および確信区間を出してくれる関数がある library(tidybayes) 点推定値 もちろんpost_sigmaに対しても同様に計算可能です EAP（事後期待値） mean(post_mu) [1] 6.98601 MED（事後中央値） median(post_mu) [1] 6.985822 MAP（事後最頻値） Mode(post_mu) [1] 6.974403 確信区間 ※乱数をもとに計算しているので 結果は微妙に異なると思います [,1] 95% HDI hdi(post_mu) 90% ETI qi(post_mu, .width = 0.9) [,2] [1,] 6.72915 7.240938 [,1] 下限と上限 [,2] [1,] 6.76985 7.198134 quantile interval とも呼ばれるのでqi 05 基本的なベイズ推論(2) 96
完成したstanコード data { int N; array[N] real SALES; } parameters { } model { どんな形のデータ(𝑌)が与えられるかを指定する。今回の例では • 100店舗の利益 (dat$sales) が与えられている。 今後データが増えたときのために 意味がわかるような名前にしておきます 推定するパラメータ(𝜃)を指定する。今回の例では • 平均 （𝜇） • 標準偏差（𝜎） の2つを推定する。 実際に事後分布の形を規定するもの(𝑃 𝑌 𝜃 𝑃 𝜃 )を指定する。 すなわち事前分布𝑃 𝜃 と尤度𝑃 𝑌 𝜃 の両方を書くブロック。 𝑃 𝜃𝑌 = } 𝑃 𝑌𝜃 𝑃 𝜃 𝑃 𝑌 05 基本的なベイズ推論(2) 97
完成したstanコード data { int N; array[N] real SALES; } parameters { } model { どんな形のデータ(𝑌)が与えられるかを指定する。今回の例では • 100店舗の利益 (dat$sales) が与えられている。 推定するパラメータ(𝜃)を指定する。今回の例では • 平均 （𝜇） • 標準偏差（𝜎） の2つを推定する。 実際に事後分布の形を規定するもの(𝑃 𝑌 𝜃 𝑃 𝜃 )を指定する。 すなわち事前分布𝑃 𝜃 と尤度𝑃 𝑌 𝜃 の両方を書くブロック。 𝑃 𝜃𝑌 = } 𝑃 𝑌𝜃 𝑃 𝜃 𝑃 𝑌 05 基本的なベイズ推論(2) 98
parametersブロック 今回推定パラメータは２つ パラメータが複数ある場合でも順番に書けば良いだけ parameters { A A real real <lower=0> mu; 平均パラメータ 𝜇 はすべての実数を取りうる sigma; 標準偏差パラメータ 𝜎 は0以上の実数を取りうる } 05 基本的なベイズ推論(2) 99
完成したstanコード data { int N; array[N] real SALES; } parameters { real mu; real <lower=0> sigma; } model { どんな形のデータ(𝑌)が与えられるかを指定する。今回の例では • 100店舗の利益 (dat$sales) が与えられている。 推定するパラメータ(𝜃)を指定する。今回の例では • 平均 （𝜇） • 標準偏差（𝜎） の2つを推定する。 実際に事後分布の形を規定するもの(𝑃 𝑌 𝜃 𝑃 𝜃 )を指定する。 すなわち事前分布𝑃 𝜃 と尤度𝑃 𝑌 𝜃 の両方を書くブロック。 𝑃 𝜃𝑌 = } 𝑃 𝑌𝜃 𝑃 𝜃 𝑃 𝑌 05 基本的なベイズ推論(2) 100
modelブロック 事前分布どうする？ 先程まで見てきた共役事前分布の例を素直に書けば p. 87 (𝑛0 = 𝜈0 = 0.001, 𝜇0 = 0, 𝜎0 = 1) model { sigma ~ inv_gamma(0.0005, 0.0005); mu ~ normal(0, sigma/sqrt(0.001)); SALES ~ normal(mu, sigma); // 尤度 } 𝑃 𝜇, 𝜎 = 𝑃 𝜇 𝜎 𝑃 𝜎 2 𝜎 𝜈0 𝜈0 𝜎02 = 𝑁 𝜇0 , 𝐼𝑛𝑣𝐺𝑎𝑚𝑚𝑎 , 𝑛0 2 2 stanの場合共役であるかどうかはあまり気にしない というかこの場合逆ガンマ分布を置くことは非推奨のレベル pp. 04 32-33 資料 資料04 pp. 32-33 共役事前分布じゃなくていいならもっと自由に，都合の良い分布を選ぼう 05 基本的なベイズ推論(2) 101
コーシー分布｜Cauchy distribution 正負の実数値を取るパラメータの事前分布として用いられる 同じ用途で𝑡分布も用いられる事が多い 自由度１の𝑡分布はコーシー分布と同じ 期待値と分散が存在しない不思議な分布 ▼実際に乱数を出してみましょう 正規分布と同じ対称形だが，裾が重い max(rcauchy(100000,0,1)) 稀にとんでもない値を出す事がある 𝑁(0,1) 𝐶𝑎𝑢𝑐ℎ𝑦(0,1) ←のような分布からたまに-1000とか出てもおかしくない 外れ値のあるモデルや回帰係数の 事前分布としても用いられる 「普通に考えるとこんなもんだけど， もしかしたらとんでもない値も無きにしもあらず」 05 基本的なベイズ推論(2) 102
事前分布を置く sigmaの事前分布 半コーシー分布をおいてみます 𝐶𝑎𝑢𝑐ℎ𝑦(0,1) 𝐶𝑎𝑢𝑐ℎ𝑦(3,1) 𝐶𝑎𝑢𝑐ℎ𝑦(0,2) sigma ~ cauchy(0, 10); • 2つのパラメータは基本的に 正規分布と同じような解釈でOK • parametersブロックで値域に制約 をかけておけば自動的に 半コーシー分布として扱われる 正則にするために 確率密度が調整されている 𝐻𝑎𝑙𝑓 − 𝐶𝑎𝑢𝑐ℎ𝑦(0,1) muの事前分布 𝐶𝑎𝑢𝑐ℎ𝑦(0,1) ふつうの正規分布でOK mu ~ normal(0, 100); 弱情報 05 基本的なベイズ推論(2) 103
完成したstanコード data { int N; array[N] real SALES; } parameters { real mu; real <lower=0> sigma; } model { mu ~ normal(0, 100); どんな形のデータ(𝑌)が与えられるかを指定する。今回の例では • 100店舗の利益 (dat$sales) が与えられている。 推定するパラメータ(𝜃)を指定する。今回の例では • 平均 （𝜇） • 標準偏差（𝜎） の2つを推定する。 実際に事後分布の形を規定するもの(𝑃 𝑌 𝜃 𝑃 𝜃 )を指定する。 すなわち事前分布𝑃 𝜃 と尤度𝑃 𝑌 𝜃 の両方を書くブロック。 sigma ~ cauchy(0, 10); SALES ~ normal(mu, sigma); } 𝑃 𝜃𝑌 = 𝑃 𝑌𝜃 𝑃 𝜃 𝑃 𝑌 05 基本的なベイズ推論(2) 104
完成したモデル
data {
model_normal.stan
int N;
array[N] real SALES;
}
parameters {
real mu;
real <lower=0> sigma;
}
model {
mu ~ normal(0, 100);
sigma ~ cauchy(0, 10);
推定に必要な情報
データ
𝑌
推定したい
パラメータ
𝜃
尤度
𝑃(𝑌|𝜃)
事前分布
今回の事例
各店舗の利益
𝑥 = (6.52,8.53, ⋯ )
利益の平均
𝜇
利益の分散
𝜎2
正規分布
𝑁𝑜𝑟𝑚𝑎𝑙(𝑥|𝜇, 𝜎 2 )
𝜇
正規分布
𝜎2
半コーシー分布
𝑃(𝜃)
【今回のベイズモデリングのplate notation】
パラメータは
2つとも
全データで共通
SALES ~ normal(mu, sigma);
𝜇
𝑥𝑖
𝜎
𝑖 data
}
05 基本的なベイズ推論(2)
105
完成したモデル
data {
model_normal.stan
int N;
array[N] real SALES;
}
parameters {
real mu;
推定に必要な情報
データ
𝑌
推定したい
パラメータ
𝜃
尤度
𝑃(𝑌|𝜃)
事前分布
𝑃(𝜃)
今回の事例
各店舗の利益
𝑥 = (6.52,8.53, ⋯ )
利益の平均
𝜇
利益の分散
𝜎2
正規分布
𝑁𝑜𝑟𝑚𝑎𝑙(𝑥|𝜇, 𝜎 2 )
𝜇
正規分布
𝜎2
半コーシー分布
real <lower=0> sigma;
}
model {
library(cmdstanr)
mu ~ normal(0, 100);
model <- cmdstan_model("model_normal.stan")
sigma ~ cauchy(0, 10);
stan_data <- list(N=100, SALES=dat$sales)
SALES ~ normal(mu, sigma);
result <- model$sample(data = stan_data)
}
05 基本的なベイズ推論(2)
106
結果を見ていく 【事後分布の要約】 その他の点推定・区間推定の方法は result$summary() 資料04で確認してください パラメータが 複数あっても 大丈夫 左から 事後平均値(EAP)，事後中央値（MED），事後分布のSD 90%確信区間(Equal-tailed interval) 【事後分布をまとめてプロット】 𝜎 2 ではなく 𝜎 の事後分布 result$draws() |> mcmc_dens(pars = c("mu", "sigma")) 05 基本的なベイズ推論(2) 107
まとめと次回予告 【まとめ】 ポアソン分布・正規分布におけるベイズ推定を行いました 共役事前分布を使えば基本的な流れは同じです ただ，正規分布レベルでも事後分布の導出は結構たいへんです stanで行う場合はそこまで気にすることはありません 【次回予告】 stanの中身（MCMC）の考え方の話をはさみます 知らなくてもベイズ推定自体はできますが， 知っていると結果のチェックやトラブルシューティングの役に立つと思います。 05 基本的なベイズ推論(2) 108