ベイズ統計_01_イントロダクション

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

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

ベイズ統計はどんどん増加している 【ざっくり検索】 Web of Scienceで • タイトル・アブストラクト・キーワードに “Bayes(ian)”を含む論文 • 分野等の指定はなし（全分野） 01 イントロダクション 3

ベイズ統計とは何なのか 01 イントロダクション 10

Q. そもそも統計学は何をしたいのか A. 不確実性を扱いたい 例 日本人の平均身長が知りたい 【学部までに習ったであろう統計学の場合】 ベイズ統計に対して「頻度主義統計学」と呼ばれたりします 理想 日本人全員に身長を聞いて，平均値を計算したら良い 現実 そんなことは無理なので，ランダムサンプリングを行って標本平均から推定する ここに不確実性がある 100人に聞いた結果，平均身長は168.2cmでした。 よって，日本人の平均身長は168.2cmです。 今回はたまたま168.2cmだったけど， 別の100人をサンプリングしたら結果変わるよね？ 01 イントロダクション 11

不確実性に基づく統計的推測 ▌統計的推測では，母集団分布の性質を調べたい 母集団分布 身長の分布 我々にはわからないが （神の視点があるとしたら） 決まった一つの分布が存在している ▶ この分布の平均値が知りたい 私達が実際に 分かる範囲 ここからどうやって 母集団分布の性質を推定する？ 標本平均 𝑥ҧ でも手元には たまたま得られた標本に基づく 統計量しかない 01 イントロダクション 12

不確実性を踏まえて推測するために 2 1 母集団の分布を なにか仮定する 出現しうる標本の パターンが分かる 3 各標本で平均値 などを計算する 標本の平均値 170.4 母集団分布 169.3 population distribution 4 標本の平均値の 分布を作れる ︙ ︙ ︙ 𝑁 𝜎2 𝜇, 𝑛 ︙ 171.4 𝑁(𝜇, 𝜎 2 ) 172.3 標本分布 sampling distribution 無限の母集団から 無限回サンプリングすると考える ここでのポイントは 無限回の試行における頻度の極限を「確率」として 標本統計量の確率分布を求めたこと 01 イントロダクション 13

不確実性に基づく統計的推測（点推定） ▌標本分布を仲介することで 母集団分布が 𝑁(𝜇, 𝜎 2 ) の場合 標本分布は 𝑁 𝜎2 𝜇, 𝑛 になる 「中心極限定理」覚えていますか？ 私達が実際に 分かる範囲 標本分布 母集団分布 標本平均 𝑥ҧ 𝜎2 𝑥,ҧ 𝑛 2 標本分布𝑁 を生み出す母集団分布は 𝑁(𝑥,ҧ 𝜎 )と考えるのが最も妥当 ▶ 母平均は 𝑥ҧ と考えるのが妥当だろう！ 01 イントロダクション 標本平均 𝑥ҧ を生み落とした 標本分布は𝑁 𝜇 = 𝜎2 𝑥,ҧ 𝑛 が最もしっくり来る 𝜎 2 がすでに分かっているとしたら 14

不確実性に基づく統計的推測（区間推定） グレーの分布＝真の母集団分布 𝑁(𝜇, 𝜎 2 ) における標本平均の標本分布 母平均の値が何であろうと 1. 標本平均の値を母平均の代わりに 使って標本分布をつくる 2. その標本分布において 95%区間を算出する を全ての標本で行うと，その区間は 必ず95%の割合で母平均を含む 全標本中5%が 作る区間は 母平均を含まない この区間を confidence interval (CI) 95%信頼区間 と呼びました 01 イントロダクション 15

別の視点から「確率」について考えてみましょう 例 「明日の降水確率は70%です。」という文の意味は？ ▌頻度主義に基づいて（無理やり）考えようとすると 母集団の分布を なにか仮定する 母集団分布 population distribution 2 出現しうる標本の パターンが分かる 標本での 生確率の 分布を作れる 晴れ ︙ 生確率 4 生確率 1 れ れ 雨 標本分布 sampling distribution 無限の母集団から 無限回サンプリングすると考える 01 イントロダクション 16

わかりにくいですね 別の問題を考えてみましょう 例 「明日の降水確率は です。」 ▌そもそも「無限回」は考えにくい そのせいで信頼区間の意味も 頻度主義統計学に基づいて （無理やり）考えようとすると 統計的仮説検定の 𝑝 値もわかりにくくなっているのです 出現しうる標本の パターンが分かる 母集団の分布を なにか仮定する 標本での 生確率の 分布を作れる 生確率 生確率 晴れ 雨 無限の母集団から 無限回サンプリングすると考える 例 データの性質 「明日の降水確率は70%です。」 つまり無限個のパラレルワールドを 考えたとしたら，という感じ 頻度主義的な解釈 同じような気圧配置が 無限回出現したときに そのうち70%では 雨が降ると思われる 01 イントロダクション 一般人の解釈 70%っていうくらいだから きっと明日は雨が 降るんでしょうな 17

一般人の思う「確率」は解釈が違う ▌我々は確率を「信念」として捉えがち 予報 頻度主義的な解釈 主観的な解釈 「明日の降水確率は10%です」 同じような条件が無限回繰り返されたときに， そのうち10%では雨が降る。 ただし明日降るかどうかは決定事項だ。 たぶん降らない。 「明日の降水確率は50%です」 同じような条件が無限回繰り返されたときに， 降るかもしれないし そのうち50%では雨が降る。 降らないかもしれない。 ただし明日降るかどうかは決定事項だ。 「明日の降水確率は90%です」 同じような条件が無限回繰り返されたときに， そのうち90%では雨が降る。 ただし明日降るかどうかは決定事項だ。 01 イントロダクション たぶん降る。 18

人間は確率を主観的なものとして捉えがち ▌頭の中には確率分布が？ 予報 「明日の降水確率は10%です」 主観的な解釈 確率分布 たぶん降らない。 れ 「明日の降水確率は50%です」 降るかもしれないし 降らないかもしれない。 れ 「明日の降水確率は90%です」 たぶん降る。 一般人の（主観的な） 解釈では 確率は「信念」のようなもの と言える気がする ベイズ統計では確率を 主観的な「信念」 として捉えている れ 01 イントロダクション 19

主観的な信念は更新される ▌情報を与えることで確率が更新される 予報 翌日の空 主観的な解釈 結構明るい 降らないかも？ 確率分布 50% 確率分布 れ 次の日 くもり どうなんでしょ。 れ れ 今にも降り出しそうな 暗い雲 やっぱり降りそう。 れ 01 イントロダクション 20

ベイズ統計学 ▌確率を主観的なものとしてあつかうと考えたときの理論 ベイズの定理 𝑃 𝑋 𝜃 𝑃(𝜃) 𝑃 𝜃𝑋 = 𝑃(𝑋) 尤度 × 事前確率 事後確率 = 周辺尤度 前ページの例に当てはめると… 𝜃：「雨が降る」 𝑃(𝜃)：天気予報による降水確率 𝑋：翌日の空模様 𝑃 𝑋 𝜃 ：雨が降るときの空模様の出現確率 雨が降る日の空は朝から暗いことが多い 雨が降らない日は空が明るいことが多い などなど… 𝑃 𝜃 𝑋 ：翌日の空模様を見た上で「雨が降る」とどの程度思っているか 事前確率をデータ（尤度）によって更新することで 不確実性を減らしていく 01 イントロダクション 21

平均身長をベイズ的に推測すると どうせ「正解」はわからないので 分布の平均値を確率分布として 考えていく 身長の 母平均に対して もともと持っている 信念 𝑃(𝜇) 手元にある情報をもとに 信念を更新していく 私達が実際に 分かる範囲 標本の値 𝒙 ▶𝑃 𝑋 𝜇 事前分布 大体170くらいだと思うけど… いくら雑に見積もっても 140以下とか200以上のはずはないよね 100人に聞いた結果， 平均身長は168.2cm，SDは5.6cmでした。 01 イントロダクション データによって 更新された 信念 𝑃 𝜇𝑋 事後分布 得られた標本も含めて考えると 母平均は大体166から170くらいの 間にあるだろう 22

なぜベイズ統計が注目されているのか 理由はいろいろあるのですが… 01 イントロダクション 23

理由１：確率の解釈が人類にあっている ▌そもそも人類は確率を頻度主義的には考えていない 信頼区間ではなく「確信区間」などと呼ばれます 【例：区間推定】 頻度主義の場合 ベイズ統計の場合 「𝜇が区間内にある確率」的な解釈は不可能 「𝜇が区間内にあるという確信の強さ」 的な解釈が可能 同じ計算を無限個の標本で繰り返す場合を考えている ▶ 手元の単一の標本に対しては何を意味している？ 01 イントロダクション 得られた事後分布自体が （標本によって更新された）信念のため 24

理由２：統計的仮説検定の問題をクリアできそう ▌統計的仮説検定の非対称性 【例：２群の差の母平均 𝜇𝐴 − 𝜇𝐵 = 𝜇𝑑 に関する検定】 頻度主義の場合 もちろんいずれにせよ 設定したモデルが正しいという 仮定のもとでの話ですが… ベイズ統計の場合 帰無仮説 H0 : 𝜇𝑑 = 0 が正しい という仮定のもとで標本を見る 𝑝 > 0.05 の場合 帰無仮説は棄却されない 𝜇𝑑 に関する信念が 確率分布として得られる 帰無仮説の正しさについては 何も言及できない H0 : 𝜇𝑑 = 0の確信度の強さを 直接計算できる 01 イントロダクション 25

（補足）なぜ「帰無仮説が正しい」といえないのか？ ▌棄却できない帰無仮説の値は１つではない 例えば帰無仮説 H0 : 𝜇𝑑 = 0 のもとで 𝑝 > 0.05であったとすると， きっと帰無仮説 H0 : 𝜇𝑑 = 0.01 のもとでも 𝑝 > 0.05 になる 厳密には，標本から作った95%信頼区間に含まれる値を帰無仮説として設定した場合は， すべて𝑝 > 0.05 となります（両側検定の場合） ▶ 「帰無仮説が棄却されない」ことを「帰無仮説が正しい」と表現すると 「𝜇𝑑 は0かつ0.01かつ0.02かつ…」となってしまう 頻度主義では母数は定数として存在するはずなので この考え方はおかしい 01 イントロダクション 26

理由３：手持ちの情報が使える 例 そんなデータを使うな，というのが 最も正しい反応ですが，あくまでも一例として… 日本人の平均身長を推定しようと思っていたのですが， 時間が無くバスケ部の10人にしか聞けませんでした。標本平均は185.2cmでした。 頻度主義の場合 ベイズ統計の場合 データから推定するしかない ◀ 最尤推定値は 𝜇Ƹ = 185.2 なので過大推定すぎる 事前確率で推定を補正することもできる 平均身長なんて きっと170くらいだろう 実際には過去のデータなど もっと妥当なやり方で事前確率を決めます 事前確率 𝑃 𝜇 = 𝑁 170, 22 事後確率 𝑃 𝜇|𝑋 = 𝑁 175.8, 1.572 いくらなんでも 平均身長が180以上 なんてことは無いだろう データによって更新 ▌もちろん，事前情報を使わないことで客観性を重視することも可能です。 01 イントロダクション 27

理由４：事後分布を数値的に求める手法が確立した ▌ 複雑な問題では事後分布を解析的に求めるのはかなり大変 ベイズの定理 𝑃 𝑋 𝜃 𝑃(𝜃) 𝑃 𝜃𝑋 = 𝑃(𝑋) （詳細は省略しますが） 右辺の分母には多重積分が登場する ▶ 実際には限られた（簡単な）場面でしか使えない理論と考えられていた ▌ 画期的な手法が研究されてきた（1990年代ごろから？） それがMCMC（マルコフ連鎖モンテカルロ法） 簡単に言えば「結果的に直接 𝑃 𝜃 𝑋 から乱数を生成したことにする」という方法 ▶ 乱数をいっぱい集めたら事後分布𝑃 𝜃 𝑋 を簡単に再生できる！ ▌ そしてMCMCを簡単に実行するソフトウェアも普及した ▶ その１つがstan ▌ その結果，複雑な統計モデリングでも割と自由に扱えるようになりました。 頻度主義の場合，標本分布を解析的に導出したりするのは大変なのです。 01 イントロダクション 28

ということで ▌ベイズはどんどん広がっているのです 豊田(2015)「基礎からのベイズ統計学」 信じるか信じないかはあなた次第です 01 イントロダクション 29

stanとは何か 余裕があればインストールまで 01 イントロダクション 30

stanって何？ ▌MCMCによるサンプリングをするための言語です 確率的プログラミング言語と呼ばれるものです。 ▌開 が盛んなので，今ひとつ始めるならこれです 他に有名なソフトウェアとしては(Multi)BUGS，JAGS, PyMCあたりでしょうか (Strumbelj et al., 2024) ▌様々な言語から呼び出すことができます R, Python, Julia, Matlab, Stata, Mathematica, Scala, そしてコマンドラインからも ▌C++で実行するので速いです stanの記法で書けば，自動的にC++に翻訳してくれます ▌あとは使えばなんとなくわかります，たぶん 最近急に赤から青に変わりました 01 イントロダクション 31

• https://www.multibugs.org/
• https://mcmc-jags.sourceforge.io/
• https://www.pymc.io/welcome.html
• https://doi.org/10.1214/23-STS907

R(またはPython)でstanを使う場合 ▌大きく分けて２種類の方法があります Rstan (Pystan) 本講義では，基本的にこちら CmdstanR (CmdstanPy) ▶ R(Python)の内部でstanを動かす ▶ R(Python)の外側で(cmd)stanを動かす 基本関数 他のパッケージ 基本関数 他のパッケージ 【メリット】 事後処理などが扱いやすい 【デメリット】 内部との干渉がいろいろあるらしい ▶ stanのバージョン更新が遅い 【メリット】 stanのバージョン更新が速い 動作も比較的安定している 【デメリット】 事後処理（結果の保存など）が少し手間 裏を返せば，そのあたりもより柔軟に設定可能です 01 イントロダクション 32

早速インストールしてみましょう (https://mc-stan.org/cmdstanr/articles/cmdstanr.html) 1. cmdstanrのインストール install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos"))) 2. C++ コンパイラのインストール 3. Cmdstanのインストール 4. 完了！ 01 イントロダクション 33

• https://mc-stan.org/cmdstanr/articles/cmdstanr.html

早速インストールしてみましょう (https://mc-stan.org/cmdstanr/articles/cmdstanr.html) 1. cmdstanrのインストール 2. C++ コンパイラのインストール いくつかやり方はありますが，Rの場合RToolsをいれるのが一般的です このページからRのバージョンにあったものRToolsをダウンロード＆インストール ここのチェックは そのまま入れておいたほうが良い 3. Cmdstanのインストール 4. 完了！ 01 イントロダクション 34

• https://cran.r-project.org/bin/windows/Rtools/
• https://mc-stan.org/cmdstanr/articles/cmdstanr.html

早速インストールしてみましょう (https://mc-stan.org/cmdstanr/articles/cmdstanr.html) 1. cmdstanrのインストール 2. C++ コンパイラのインストール 3. Cmdstanのインストール まずはインストールの準備が整っているか確認 library(cmdstanr) check_cmdstan_toolchain() The C++ toolchain required for CmdStan is setup properly! となればOK！ 4. 完了！ 01 イントロダクション 35

• https://mc-stan.org/cmdstanr/articles/cmdstanr.html

どうやらC++コンパイラのインストールがうまく行かない場合は… ▌とりあえずRのバージョンを最新にしてください。 ▌（たぶん）バージョンがある程度新しければ check_cmdstan_toolchain() 後に Run cmdstanr::check_cmdstan_toolchain(fix = TRUE) to fix the issue. 的なことが表示されます。 ▌表示されたらそれに従って check_cmdstan_toolchain(fix = TRUE) を実行しましょう。 ▶ 詳細はわからないのですが，とりあえずcmdstanのインストールに必要なものを うまいこと入れてくれるようです。 01 イントロダクション 36

早速インストールしてみましょう (https://mc-stan.org/cmdstanr/articles/cmdstanr.html) 1. cmdstanrのインストール 2. C++ コンパイラのインストール 3. Cmdstanのインストール いよいよインストールを実行 install_cmdstan(cores = 4) # coresは自分のPCのコア数に合わせて変えてください （そこそこ時間がかかると思います。） 4. 完了！ 01 イントロダクション 37

• https://mc-stan.org/cmdstanr/articles/cmdstanr.html

早速インストールしてみましょう (https://mc-stan.org/cmdstanr/articles/cmdstanr.html) 1. cmdstanrのインストール 2. C++ コンパイラのインストール 3. Cmdstanのインストール 4. 完了！ 正しくcmdstanが呼び出せるか確認 cmdstan_path() "/Users/***/.cmdstan/cmdstan-***" 01 イントロダクション と表示されればたぶんOK！ 38

• https://mc-stan.org/cmdstanr/articles/cmdstanr.html

どうやらcmdstanのインストールがうまく行かない場合は… cmdstanをインストールしようとしているディレクトリの問題かもしれません。 ▌cmdstanインストール時に，cmdstanの保存先を指定しましょう。 （一例） install_cmdstan(dir = "C:/Users/Bunji/hoge/fuga”, cores = 4) （一案） install_cmdstan(dir = "C:/cmdstan”, cores = 4) 自分のPCの中の cmdstanが入っていても じゃまにならない場所 ここで指定するフォルダは 既存のフォルダ，または事前に「新しいフォルダ」で 作っておく必要があります。 うまく行けば，指定したフォルダ内に ”cmdstan-x.yy.z”（後ろはバージョン番号）ができるはず そして 一応確認 01 イントロダクション 39

インストールが完了したら ▌うまくいくかを確認しておいてください 以下のコードをまるまるコピーして実行してみましょう file <- file.path(cmdstan_path(), "examples", "bernoulli", "bernoulli.stan") mod <- cmdstan_model(file) data_list <- list(N = 10, y = c(0,1,0,0,0,0,0,0,0,1)) fit <- mod$sample(data = data_list, chains = 2) （前略） Chain 1 Iteration: 1900 / 2000 [ 95%] (Sampling) Chain 1 Iteration: 2000 / 2000 [100%] (Sampling) Chain 1 finished in 0.0 seconds. 01 イントロダクション という感じで （エラー等なく） 表示されれば ほぼ確実にOK！ 40

どうやらRパッケージのインストールがうまく行かない場合は… ▌Windowsの場合，ディレクトリの問題かもしれません。 ▌最も多いのは「ドキュメント」フォルダがOneDriveと同期されている場合 install.packages()でインストールされたRのパッケージは，デフォルトでは「ドキュメント」フォルダ内 （例）C:¥Users¥Bunji¥Documents¥R¥win-library¥4.x に保存されるのですが，何も考えずに（デフォルト設定で）OneDriveを同期させると，このドキュメントフォルダが勝手に （例）C:¥Users¥Bunji¥OneDrive – 神戸大学【全学】 みたいな名前に置き換わってしまいます。Rはまだ日本語に強くないので，このようにフォルダ名に日本語（2バイト）文字が入っ ているとうまくいかないのです。 ▶ 修正はこのページを参照すると良いと思います。 （ただしレジストリを弄る必要があるのでやや難しいかも） ▌「パソコンが苦手な人」のために，次ページでは別の方法を紹介します。 Rを別の場所にインストールする試みです。 01 イントロダクション 41

• https://norimune.net/3609

Rのインストール先を指定する Rのインストールを改めてやり直す過程で… ▌「インストール先の指定」でOneDriveの管轄外のフォルダを指定する （例）Program Filesフォルダはたぶんいいぞ 他にも設定が必要かもしれないので， うまく行かなければご連絡ください。 あとはいつもと同じように インストールしてみてください。 忘れずに！ 特にこだわりがなければ C:¥Program Files¥R¥R-4.x.x は誰もが問題なく使えるはず 01 イントロダクション ディレクトリの一番うしろ（バージョン番号）は インストールしようとしているRのバージョンと 揃えておいてください 42