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

16.5K Views

April 09, 24

スライド概要

神戸大学大学院経営学研究科で2024年度より開講している「ベイズ統計」の講義資料「01_イントロダクション」を公開用に調整したものです。ベイズ統計とは何かを簡潔に紹介し,なぜベイズ統計が利用されているのか,その理由をいくつか紹介しています。また,cmdstanrのインストールも簡単に解説しています。

profile-image

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

シェア

またはPlayer版

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

関連スライド

各ページのテキスト
1.

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

2.

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

3.

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

4.

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

5.

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

6.

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

7.

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

8.

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

9.

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

10.

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

11.

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

12.

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

13.

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

14.

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

15.

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

16.

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

17.

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

18.

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

19.

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

20.

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

21.

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

22.

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

23.

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

24.

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

25.

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

26.

早速インストールしてみましょう (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

27.

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

28.

早速インストールしてみましょう (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

29.

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

30.

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

31.

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

32.

どうやら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

33.

インストールが完了したら ▌うまくいくかを確認しておいてください 以下のコードをまるまるコピーして実行してみましょう 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

34.

どうやら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

35.

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