医療職からデータサイエンティストへ医療職からデータサイエンティストへ統計学、機械学習に関する記事をまとめています。ベイズ今回は、Rのstanを使ってベイズ推定を使った単回帰分析を行なっていきます。本来であればベイズ推定を使わなくても単回帰分析のパラメーターは推定できるのですが、stanに慣れるためにもまずは簡単なところからですね。最後には通常の単回帰分析と結果の比較も行います。今回のデータは生徒100人の一日の勉強時間とテスト点数の仮想データとしましょう。geom_point()目的変数Yがテストの点数です。このデータに通常の回帰分析を行います。2.2e-16という結果になりました。1日の勉強時間が1時間増えるとテストの点数が7.2点上昇するようです。先ほどのグラフに回帰直線の95%信頼区間と予測値の信頼区間を書いてみます。title("lm_model")青色が回帰直線の信頼区間、破線が予測値の信頼区間を表しています。ここまでは普通の単回帰分析ですね。何をするかさて、stanを使ってベイズ推定を行う前に、ベイズ推定の流れを説明しておきます。ベイズ推定では以下のように単回帰分析を行います。従う正規分布として求める。←予測値の確率分布を求める。以上がベイズ推定の流れになります。やるかそれでは、ざっくりと大枠を掴んだところで、さっそくRとstanを使ってモデリングしていきます。1ですが、右辺の分母は定数と扱えるのでとしておきます。なる確率を表すので尤度関数になります。正規分布に従うのでと表すことができます。この正規分布から得られる確率密度の1からnまでの掛け合わせが尤度となるので、先ほどの尤度関数はとなります。パラメーターの事前分布ですが、今回は何も情報がないので無情報分布としておきましょう。すると結局をモデリングすることになります。これで1の準備が整いました。stanでモデリングしてみます。parameteresはdataやparametersで定義されている値を使って新たにサンプリングを行います。今回は、サンプリングされたa,bとデータstudyを使って回帰直線muをサンプリングしています。掛け合わせ(厳密には内部で対数尤度を足している)を行っています。本来であれば、事前分布の分布も定義する必要がありますが、stanではパラメーターの分布を何も指定しなければ、勝手に無情報分布を選んでくれます。ここまでで1ができました。続いて2をモデリングしていきましょう。"prediction")に相当することを行なっていきます。これをモデリングするためにstanファイルに以下を書き加えます。}}まず、最初のdataにstudy_newを渡すための変数を加えます。quantitiesは、サンプリング後のパラメーターを使って新たな値を作成します。parameteresと違って、全てのサンプリングが終了した後に実行されます。そのため、stanファイルに書かなくても後から作成することができるのですが、こちらで定義しておいた方がスムーズに結果を出すことができるでしょう。これで準備が整いました。最終的なstanコードは以下になります。}}stanは上から順に実行されるので記述の順番に気をつけてください。このstanファイルをregression.stanとして保存します。rstan::extract(stan_fit_res)さてこれでベイズ推定が完了しました。回帰直線の切片と傾きを通常の回帰分析と比べてみましょう。同じ結果になっていますね。続いて回帰直線の95%信頼区間、予測値の95%信頼区間をグラフにしてみましょう。扱いやすいようにデータフレームに加工します。91.85203これを先ほどと同じようにグラフ化します。どこかでみたようなグラフですね。これはパラメーターの事前分布に無情報分布を選んだため、通常の最尤法で求めた結果とほぼ同じになります。まとめ今回はベイズ推定で単回帰分析行い、その結果を通常の単回帰分析と比較しました。ベイズ推定では、事前分布に無情報分布を選ぶと通常の最尤法と同じ結果になります。前回も書いたようにベイズ推定の良いところは、さらに複雑なモデリングができるところにあります。stanにもだいぶ慣れてきたので、次回からは複雑なモデリングにも挑戦していきます!medi-data.hatenablog.com※本記事は筆者が個人的に学んだことをまとめた記事なります。数学の記法や詳細な理論等で誤りがあった際はご指摘頂けると幸いです。参考こちらの内容をアレンジして記事を書きました。これほどわかりやすくベイズモデリングを学べる書籍はほとんどないでしょう。前半にベイズ推定の流れが書いてあります。前回はベイズ推定を使って単回帰分析を行いました。学んだこと、日々の所感をまとめていきます。※本ブログは筆者が個人的に学んだこと感じたことをまとめたブログになります。所属する組織の意見・見解とは無関係です。medi_data0826引用をストックしましたストック一覧を見る閉じる引用するにはまずログインしてくださいログイン閉じる引用をストックできませんでした。再度お試しください閉じる限定公開記事のため引用できません。