35.7K Views
April 09, 24
スライド概要
神戸大学大学院経営学研究科で2024年度より開講している「ベイズ統計」の講義資料「01_イントロダクション」を公開用に調整したものです。ベイズ統計とは何かを簡潔に紹介し,なぜベイズ統計が利用されているのか,その理由をいくつか紹介しています。また,cmdstanrのインストールも簡単に解説しています。
【更新履歴】
・2025/08/04:cmdstanrのインストール情報のリンクを更新(pp. 34-39)
神戸大学経営学研究科准教授 分寺杏介(ぶんじ・きょうすけ)です。 主に心理学的な測定・教育測定に関する研究を行っています。 講義資料や学会発表のスライドを公開していきます。 ※スライドに誤りを見つけた方は,炎上させずにこっそりお伝えいただけると幸いです。
ベイズ統計 01 イントロダクション 分寺 杏介 神戸大学大学院 経営学研究科 [email protected] ※本スライドは,クリエイティブ・コモンズ 表示-非営利 4.0 国際 ライセンス(CC BY-NC 4.0)に従って利用が可能です。
授業概要・到達目標 【授業概要】 これから社会科学の研究者になるなら,ベイズ統計は多分知らないとまずいです 【到達目標】 ▌ベイズ統計学の理論的な側面について,基本的な理解を獲得する ベイズ統計を用いている論文を読むときに「何をやっているか」が分かるようになってほしいのです ▌stanを用いて,実際に自分でベイズ統計学に基づく分析ができるようになる 自分の研究で使ってみようと思ってくれる人が現れたらうれしいのです ▌ベイズ統計学における分析上の注意点を理解し,正しく運用できるようになる どうせやるなら「なんとなくできる」ではなく「きちんとできる」ようになってほしいのです 01 イントロダクション 2
ベイズ統計はどんどん増加している 【ざっくり検索】 Web of Scienceで • タイトル・アブストラクト・キーワードに “Bayes(ian)”を含む論文 • 分野等の指定はなし(全分野) ※ 01 イントロダクション 3年までのデータです 3
ベイズ統計とは何なのか 01 イントロダクション 11
Q. そもそも統計学は何をしたいのか A. 不確実性を扱いたい 例 日本人の平均身長が知りたい 【学部までに習ったであろう統計学の場合】 ベイズ統計に対して「頻度主義統計学」と呼ばれたりします 理想 日本人全員に身長を聞いて,平均値を計算したら良い 現実 そんなことは無理なので,ランダムサンプリングを行って標本平均から推定する ここに不確実性がある 100人に聞いた結果,平均身長は168.2cmでした。 よって,日本人の平均身長は168.2cmです。 今回はたまたま168.2cmだったけど, 別の100人をサンプリングしたら結果変わるよね? 01 イントロダクション 12
不確実性に基づく統計的推測 ▌統計的推測では,母集団分布の性質を調べたい 母集団分布 身長の分布 3 我々にはわからないが (神の視点があるとしたら) 決まった一つの分布が存在している この分布の平均値が知りたい 私達が実際に 分かる範囲 ここからどうやって 母集団分布の性質を推定する? 標本平均 𝑥ҧ でも手元には たまたま得られた標本に基づく 統計量しかない 01 イントロダクション 13
不確実性を踏まえて推測するために 2 出現しうる標本の 1 母集団の分布を パターンが分かる 3 各標本で平均値 などを計算する なにか仮定する 標本の平均値 母集団分布 170.4 4 標本の平均値の 分布を作れる 169.3 population distribution ︙ ︙ 3 ︙ 𝑁 3 𝜎2 𝜇, 𝑛 ︙ 171.4 𝑁(𝜇, 𝜎 2 ) 172.3 3 標本分布 sampling distribution 無限の母集団から 無限回サンプリングすると考える ここでのポイントは 無限回の試行における頻度の極限を「確率」として 標本統計量の確率分布を求めたこと 01 イントロダクション 14
不確実性に基づく統計的推測(点推定) ▌標本分布を仲介することで 母集団分布が 𝑁(𝜇, 𝜎 2 ) の場合 標本分布は 𝑁 𝜎2 𝜇, 𝑛 になる 「中心極限定理」覚えていますか? 私達が実際に 分かる範囲 標本分布 母集団分布 3 3 標本平均 𝑥ҧ 𝜎2 標本分布𝑁 𝑥,ҧ を生み出す母集団分布は 𝑛 𝑁(𝑥,ҧ 𝜎 2 )と考えるのが最も妥当 母平均は 𝑥ҧ と考えるのが妥当だろう! 01 イントロダクション 標本平均 𝑥ҧ を生み落とした 標本分布は𝑁 𝜎2 𝜇 = 𝑥,ҧ 𝑛 が最もしっくり来る 𝜎 2 がすでに分かっているとしたら 15
不確実性に基づく統計的推測(区間推定) グレーの分布=真の母集団分布 𝑁(𝜇, 𝜎 2 ) における標本平均の標本分布 母平均の値が何であろうと 1. 標本平均の値を母平均の代わりに 使って標本分布をつくる 2. その標本分布において 95%区間を算出する を全ての標本で行うと,その区間は 必ず95%の割合で母平均を含む 全標本中5%が 作る区間は 母平均を含まない この区間を confidence interval (CI) 95%信頼区間 と呼びました 01 イントロダクション 16
別の視点から「確率」について考えてみましょう 例 「明日の降水確率は70%です。」という文の意味は? ▌頻度主義に基づいて(無理やり)考えようとすると なにか仮定する 母集団分布 population distribution 2 出現しうる標本の パターンが分かる 分布を作れる 晴れ ︙ 生確率 4 標本での 生確率の 生確率 1 母集団の分布を れ れ 雨 標本分布 sampling distribution 無限の母集団から 無限回サンプリングすると考える 01 イントロダクション 17
わかりにくいですね 別の問題を考えてみましょう 例 「明日の降水確率は です。」 ▌そもそも「無限回」は考えにくい そのせいで信頼区間の意味も 頻度主義統計学に基づいて (無理やり)考えようとすると 統計的仮説検定の 𝑝 値もわかりにくくなっているのです 出現しうる標本の パターンが分かる 母集団の分布を なにか仮定する 標本での 生確率の 分布を作れる 生確率 生確率 晴れ 雨 無限の母集団から 無限回サンプリングすると考える 例 「明日の降水確率は70%です。」 データの性 つまり無限個のパラレルワールドを 考えたとしたら,という感じ 頻度主義的な解釈 同じような気圧配置が 無限回出現したときに そのうち70%では 雨が降ると思われる 01 イントロダクション 一般人の解釈 70%っていうくらいだから きっと明日は雨が 降るんでしょうな 18
一般人の思う「確率」は解釈が違う ▌我々は確率を「信念」として捉えがち 予報 頻度主義的な解釈 主観的な解釈 「明日の降水確率は10%です」 同じような条件が無限回繰り返されたときに, そのうち10%では雨が降る。 ただし明日降るかどうかは決定事項だ。 たぶん降らない。 「明日の降水確率は50%です」 同じような条件が無限回繰り返されたときに, 降るかもしれないし そのうち50%では雨が降る。 降らないかもしれない。 ただし明日降るかどうかは決定事項だ。 「明日の降水確率は90%です」 同じような条件が無限回繰り返されたときに, そのうち90%では雨が降る。 ただし明日降るかどうかは決定事項だ。 01 イントロダクション たぶん降る。 19
人間は確率を主観的なものとして捉えがち ▌頭の中には確率分布が? 予報 「明日の降水確率は10%です」 主観的な解釈 確率分布 たぶん降らない。 れ 「明日の降水確率は50%です」 降るかもしれないし 降らないかもしれない。 れ 「明日の降水確率は90%です」 たぶん降る。 一般人の(主観的な) 解釈では 確率は「信念」のようなもの と言える気がする ベイズ統計では確率を 主観的な「信念」 として捉えている れ 01 イントロダクション 20
主観的な信念は更新される ▌情報を与えることで確率が更新される 予報 翌日の空 主観的な解釈 結構明るい 降らないかも? 確率分布 50% 確率分布 れ 次の日 くもり どうなんでしょ。 れ れ 今にも降り出しそうな 暗い雲 やっぱり降りそう。 れ 01 イントロダクション 21
ベイズ統計学 ▌確率を主観的なものとしてあつかうと考えたときの理論 ベイズの定理 𝑃 𝑋 𝜃 𝑃(𝜃) 𝑃 𝜃𝑋 = 𝑃(𝑋) 尤度 × 事前確率 事後確率 = 周辺尤度 前ページの例に当てはめると… 𝜃:「雨が降る」 𝑃(𝜃):天気予報による降水確率 𝑋:翌日の空模様 𝑃 𝑋 𝜃 :雨が降るときの空模様の出現確率 雨が降る日の空は朝から暗いことが多い 雨が降らない日は空が明るいことが多い などなど… 𝑃 𝜃 𝑋 :翌日の空模様を見た上で「雨が降る」とどの程度思っているか 事前確率をデータ(尤度)によって更新することで 不確実性を減らしていく 01 イントロダクション 22
平均身長をベイズ的に推測すると どうせ「正解」はわからないので 分布の平均値を確率分布として 考えていく 手元にある情報をもとに 信念を更新していく 身長の 私達が実際に 分かる範囲 母平均に対して 3 もともと持っている 信念 𝑃(𝜇) 標本の値 𝒙 3 データによって 更新された 信念 𝑃 𝜇𝑋 𝑃 𝑋𝜇 事前分布 大体170くらいだと思うけど… いくら雑に見積もっても 140以下とか200以上のはずはないよね 100人に聞いた結果, 平均身長は168.2cm,SDは5.6cmでした。 01 イントロダクション 事後分布 得られた標本も含めて考えると 母平均は大体166から170くらいの 間にあるだろう 23
なぜベイズ統計が注目されているのか 理由はいろいろあるのですが… 01 イントロダクション 24
理由1:確率の解釈が人類にあっている ▌そもそも人類は確率を頻度主義的には考えていない 信頼区間ではなく「確信区間」などと呼ばれます 【例:区間推定】 頻度主義の場合 ベイズ統計の場合 「𝜇が区間内にある確率」的な解釈は不可能 「𝜇が区間内にあるという確信の強さ」 的な解釈が可能 3 同じ計算を無限個の標本で繰り返す場合を考えている 手元の単一の標本に対しては何を意味している? 01 イントロダクション 得られた事後分布自体が (標本によって更新された)信念のため 25
理由2:統計的仮説検定の問題をクリアできそう ▌統計的仮説検定の非対称性 【例:2群の差の母平均 𝜇𝐴 − 𝜇𝐵 = 𝜇𝑑 に関する検定】 頻度主義の場合 もちろんいずれにせよ 設定したモデルが正しいという 仮定のもとでの話ですが… ベイズ統計の場合 帰無仮説 H0 : 𝜇𝑑 = 0 が正しい という仮定のもとで標本を見る 𝑝 > 0.05 の場合 帰無仮説は棄却されない 𝜇𝑑 に関する信念が 確率分布として得られる 帰無仮説の正しさについては 何も言及できない H0 : 𝜇𝑑 = 0の確信度の強さを 直接計算できる 01 イントロダクション 26
(補足)なぜ「帰無仮説が正しい」といえないのか? ▌棄却できない帰無仮説の値は1つではない 例えば帰無仮説 H0 : 𝜇𝑑 = 0 のもとで 𝑝 > 0.05であったとすると, きっと帰無仮説 H0 : 𝜇𝑑 = 0.01 のもとでも 𝑝 > 0.05 になる 厳密には,標本から作った95%信頼区間に含まれる値を帰無仮説として設定した場合は, すべて𝑝 > 0.05 となります(両側検定の場合) 「帰無仮説が棄却されない」ことを「帰無仮説が正しい」と表現すると 「𝜇𝑑 は0かつ0.01かつ0.02かつ…」となってしまう 頻度主義では母数は定数として存在するはずなので この考え方はおかしい 01 イントロダクション 27
理由3:手持ちの情報が使える 例 日本人の平均身長を推定しようと思っていたのですが, そんなデータを使うな,というのが 最も正しい反応ですが,あくまでも一例として… 時間が無くバスケ部の10人にしか聞けませんでした。標本平均は185.2cmでした。 頻度主義の場合 データから推定するしかない ベイズ統計の場合 事前確率で推定を補正することもできる 平均身長なんて きっと170くらいだろう 最尤推定値は 𝜇Ƹ = 185.2 なので過大推定すぎる 実際には過去のデータなど もっと妥当なやり方で事前確率を決めます 事前確率 𝑃 𝜇 = 𝑁 170, 22 事後確率 𝑃 𝜇|𝑋 = 𝑁 175.8, 1.572 いくらなんでも 平均身長が180以上 なんてことは無いだろう データによって更新 ▌もちろん,事前情報を使わないことで客観性を重視することも可能です。 01 イントロダクション 28
理由4:事後分布を数値的に求める手法が確立した ▌ 複雑な問題では事後分布を解析的に求めるのはかなり大変 ベイズの定理 𝑃 𝑋 𝜃 𝑃(𝜃) 𝑃 𝜃𝑋 = 𝑃(𝑋) (詳細は省略しますが) 右辺の分母には多重積分が登場する 実際には限られた(簡単な)場面でしか使えない理論と考えられていた ▌ 画期的な手法が研究されてきた(1990年代ごろから?) それがMCMC(マルコフ連鎖モンテカルロ法) 簡単に言えば「結果的に直接 𝑃 𝜃 𝑋 から乱数を生成したことにする」という方法 乱数をいっぱい集めたら事後分布𝑃 𝜃 𝑋 を簡単に再生できる! ▌ そしてMCMCを簡単に実行するソフトウェアも普及した その1つがstan ▌ その結果,複雑な統計モデリングでも割と自由に扱えるようになりました。 頻度主義の場合,標本分布を解析的に導出したりするのは大変なのです。 01 イントロダクション 29
ということで ▌ベイズはどんどん広がっているのです 豊田(2015)「基礎からのベイズ統計学」 信じるか信じないかはあなた次第です 01 イントロダクション 30
stanとは何か 余裕があればインストールまで 01 イントロダクション 31
stanって何? ▌MCMCによるサンプリングをするための言語です 確率的プログラミング言語と呼ばれるものです。 ▌開 が盛んなので,今ひとつ始めるならこれです 他に有名なソフトウェアとしては(Multi)BUGS,JAGS, PyMCあたりでしょうか (Strumbelj et al., 2024) ▌様々な言語から呼び出すことができます R, Python, Julia, Matlab, Stata, Mathematica, Scala, そしてコマンドラインからも ▌C++で実行するので速いです stanの記法で書けば,自動的にC++に翻訳してくれます ▌あとは使えばなんとなくわかります,たぶん 最近急に赤から青に変わりました 01 イントロダクション 32
R(またはPython)でstanを使う場合 ▌大きく分けて2種類の方法があります Rstan (Pystan) 本講義では,基本的にこちら CmdstanR (CmdstanPy) R(Python)の内部でstanを動かす R(Python)の外側で(cmd)stanを動かす 基本関数 他のパッケージ 基本関数 他のパッケージ 【メリット】 事後処理などが扱いやすい 【デメリット】 内部との干渉がいろいろあるらしい stanのバージョン更新が遅い 【メリット】 stanのバージョン更新が速い 動作も比較的安定している 【デメリット】 事後処理(結果の保存など)が少し手間 裏を返せば,そのあたりもより柔軟に設定可能です 01 イントロダクション 33
早速インストールしてみましょう (https://stan-dev.r-universe.dev/articles/cmdstanr/cmdstanr.html) 1. cmdstanrのインストール install.packages("cmdstanr", repos = c("https://stan-dev.r-universe.dev", "https://cloud.r-project.org")) 2. C++ コンパイラのインストール 3. Cmdstanのインストール 4. 完了! 01 イントロダクション 34
早速インストールしてみましょう (https://stan-dev.r-universe.dev/articles/cmdstanr/cmdstanr.html) 1. cmdstanrのインストール 2. C++ コンパイラのインストール cmdstanをインストールするために必要なものがコンパイラです まずはもともとコンパイラが入っていないかチェック library(cmdstanr) check_cmdstan_toolchain() The C++ toolchain required for CmdStan is setup properly! (だめな例) となればOK! Error: Rtools** was not found but is required to run CmdStan with R version *.*.*. Please install Rtools** and run cmdstanr::check_cmdstan_toolchain(). 01 イントロダクション 35
早速インストールしてみましょう (https://stan-dev.r-universe.dev/articles/cmdstanr/cmdstanr.html) 1. cmdstanrのインストール 2. C++ コンパイラのインストール その他のOSについてはこちらを参照 コンパイラがない場合,WindowsならばRToolsを使うのが一般的です このページからRのバージョンにあったものRToolsをダウンロード&インストール ここのチェックは そのまま入れておいたほうが良い 3. Cmdstanのインストール 4. 完了! 01 イントロダクション 36
どうやらC++コンパイラのインストールがうまく行かない場合は… ▌とりあえずRのバージョンを最新にしてください。 ▌(たぶん)バージョンがある程度新しければ check_cmdstan_toolchain() 後に Run cmdstanr::check_cmdstan_toolchain(fix = TRUE) to fix the issue. 的なことが表示されます。 ▌表示されたらそれに従って check_cmdstan_toolchain(fix = TRUE) を実行しましょう。 詳細はわからないのですが,とりあえずcmdstanのインストールに必要なものを うまいこと入れてくれるようです。 01 イントロダクション 37
早速インストールしてみましょう (https://stan-dev.r-universe.dev/articles/cmdstanr/cmdstanr.html) 1. cmdstanrのインストール 2. C++ コンパイラのインストール 3. Cmdstanのインストール いよいよインストールを実行 install_cmdstan(cores = 4) # coresは自分のPCのコア数に合わせて変えてください (そこそこ時間がかかると思います。) 4. 完了! 01 イントロダクション 38
早速インストールしてみましょう (https://stan-dev.r-universe.dev/articles/cmdstanr/cmdstanr.html) 1. cmdstanrのインストール 2. C++ コンパイラのインストール 3. Cmdstanのインストール 4. 完了! 正しくcmdstanが呼び出せるか確認 cmdstan_path() "/Users/***/.cmdstan/cmdstan-***" 01 イントロダクション など,どこかしらのディレクトリが 表示されればたぶんOK! 39
どうやら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 イントロダクション 40
インストールが完了したら ▌うまくいくかを確認しておいてください 以下のコードをまるまるコピーして実行してみましょう 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! 41
どうやらRパッケージのインストールがうまく行かない場合は… ▌Windowsの場合,ディレクトリの問題かもしれません。 ▌最も多いのは「ドキュメント」フォルダがOneDriveと同期されている場合 install.packages()でインストールされたRのパッケージは,デフォルトでは「ドキュメント」フォルダ内 (例)C:¥Users¥Bunji¥Documents¥R¥win-library¥4.x に保存されるのですが,何も考えずに(デフォルト設定で)OneDriveを同期させると,このドキュメントフォルダが勝手に (例)C:¥Users¥Bunji¥OneDrive – 神戸大学【全学】 みたいな名前に置き換わってしまいます。Rはまだ日本語に強くないので,このようにフォルダ名に日本語(2バイト)文字が入っ ているとうまくいかないのです。 修正はこのページを参照すると良いと思います。 (ただしレジストリを弄る必要があるのでやや難しいかも) ▌「パソコンが苦手な人」のために,次ページでは別の方法を紹介します。 Rを別の場所にインストールする試みです。 01 イントロダクション 42
Rのインストール先を指定する Rのインストールを改めてやり直す過程で… ▌「インストール先の指定」でOneDriveの管轄外のフォルダを指定する (例)Program Filesフォルダはたぶんいいぞ 他にも設定が必要かもしれないので, うまく行かなければご連絡ください。 あとはいつもと同じように インストールしてみてください。 忘れずに! 特にこだわりがなければ C:¥Program Files¥R¥R-4.x.x は誰もが問題なく使えるはず 01 イントロダクション ディレクトリの一番うしろ(バージョン番号)は インストールしようとしているRのバージョンと 揃えておいてください 43