0$の値を指定します。事前分布を計算します。0$を作成して、値ごとに確率密度を計算します。この例では、0からlambda_truthの2倍を範囲とします。ガンマ分布の確率密度は、dgamma()で計算できます。確率変数の引数xにlambda_vec、形状の引数shapeにa、尺度の引数rateにbを指定します。事前分布のグラフを作成します。"density")事前分布(ガンマ分布)のグラフ真のパラメータの値を赤色の破線で示します。観測データ$\mathbf{X}$からパラメータ$\lambda$の事後分布を求めます(パラメータ$\lambda$を分布推定します)。\hat{b}$を計算します。101事後分布のパラメータは、次の式で計算できます。x_n$とデータ数$N$を加えています。事後分布を計算します。b_hatを使って、事前分布のときと同様にして処理します。事後分布のグラフを作成します。"density")事後分布(ガンマ分布)のグラフ$\lambda$の真の値付近をピークとする分布になっています。\hat{b}$から未観測のデータ$x_{*}$の予測分布を求めます。b$を用いて、予測分布(負の二項分布)のパラメータを計算します。0.009803922予測分布のパラメータは、次の式で計算できます。b_hatを使って計算できます。bを使って計算できます。予測分布を計算します。0.000782負の二項分布の確率は、dnbinom()で計算できます。失敗回数の引数xにx_vec、成功回数の引数sizeにr_hat、成功確率の引数probに1-p_hatを指定します。予測分布のグラフを真の分布と重ねて作成します。"probability")予測分布(負の二項分布)のグラフ観測データが増えると、予測分布が真の分布に近付きます。以上で、ポアソン分布のベイズ推論を実装できました。学習推移の可視化gganimateパッケージを利用して、パラメータの推定値(更新値)の推移(分布の変化)をアニメーション(gif画像)で確認します。lambda_truth*3)実装時と同じく、パラメータを指定して、確率変数がとり得る値(x軸の値)を作成します。x_vecの要素数(値の間隔)を調整してください。事後分布と予測分布のパラメータの計算(更新)処理について、2通りの方法を載せます。目的に応じて使い分けてください。推論処理:for関数による処理1つ目の処理方法では、for()を使って、1データずつ生成してパラメータの更新を繰り返し処理します。こちらの方が、前ステップで求めた事後分布(のパラメータ)を次ステップの事前分布(のパラメータ)として用いる逐次学習をイメージしやすいです。・コード(クリックで展開)超パラメータの初期値を使って、事前分布を計算します。rows現在のパラメータを文字列結合し因子型に変換して、フレーム切替用のラベル列とします。事前分布のパラメータを使って、予測分布のパラメータと予測分布を計算します。p=0.5同様に、フレーム切替用のラベル列を作成します。パラメータの更新処理をN回繰り返します。bを繰り返し更新(上書き)していきます。x_n$と$N$の計算は、ループ処理によってN回繰り返しx_n[n]と1を加えることで行います。n回目のループ処理のときには、n-1回分のx_n[n]と1が既にaとbに加えられています。更新したパラメータを使って事後分布と予測分布を計算して、それぞれ計算結果をanime_***_dfに結合していきます。x_vecの要素数」掛ける「N+1」行のデータフレームになります。行数が増えすぎないように注意してくださいしてください。推論処理:tidyverseパッケージによる処理2つ目の処理方法では、tidyverseパッケージの関数を使って、一度の処理でN+1回分の計算をします。こちらの方が、処理時間が短いです。・コード(クリックで展開)ポアソン分布に従う$N$個のデータを生成します。6ラベル用のテキストとして、観測データを集計します。0"試行ごとに、それまでの試行における各値の度数を求めて、文字列結合したベクトルを作成します。(0回目に対応する欠損値を含めた)観測データx_n、試行回数(データ番号)0:N、度数計算用の値1をデータフレームに格納します。$x$の値x_vecとデータ番号の全ての組み合わせをexpand_grid()で作成して、そこに観測データをright_join()で結合します。x_nにない値はfreq列が欠損値NAになるので、replace_na()で0に置換します。x列($x$の値)でグループ化して、freq列(各試行における観測値が1でそれ以外は0)の累積和をcumsum()で計算して、各試行回数までの度数を求めます。pivot_longer()で、試行回数列と、x_vecに対応した度数列に展開します。度数列の値をunite()で文字列結合して、pull()でベクトルとして取り出します。事前分布を含めたN+1回分の事後分布を計算します。rows試行回数(観測データ数)を表す0からNまでのN+1個の整数と、x軸の値lambda_vecの全ての組み合わせをexpand_grid()で作成します。これにより、lambda_vecの各要素をN+1フレーム分に複製できます。事前分布とN回分の事後分布のパラメータa(のベクトル)を計算(作成)して、n列の値をインデックスとして使って各試行に対応する値を抽出します。bについてはn列を使って計算できます。確率変数とパラメータの組み合わせごとに確率密度を計算して、パラメータごとにフレーム切替用のラベルを作成します。因子型への変換処理では、無名関数function()の省略記法\()を使って、(\(x){引数xを使った具体的な処理})()としています。直前のパイプ演算子を%>%にすると、行全体(\(引数){処理})()を{}の中の処理に置き換えられます(置き換えられるように引数名を.にしています)。同様に、初期値を含めたN+1回分の予測分布を計算します。rows0からNまでの整数と、x軸の値x_vecの全ての組み合わせを作成します。N+1回分の予測分布パラメータr(のベクトル)を計算(作成)して、各試行に対応する値を抽出します。pについてはn列を使って計算できます。予測分布の推移をそれぞれアニメーションで可視化します。・コード(クリックで展開)観測データと対応するラベルをデータフレームに格納します。rows事前分布には観測データが影響しないので欠損値NAとして、anime_posterior_dfのラベル列と合わせて格納します。参考として、各試行における観測データを表示するのに使います。事後分布の推移のアニメーションを作成します。600)フレームの順番を示す列をtransition_manual()に指定して、animate()でgif画像を作成します。事前分布(初期値)を含むため、フレーム数はN+1です。真の分布を複製して、対応するラベルとデータフレームに格納します。フレーム切替用ラベル)anime_model_df棒グラフの場合は、フレームごとにデータを用意する必要があります。真の分布は全てのフレームで変化しないので、N+1個に複製します。観測データと対応するラベルをデータフレームに格納します。rowsこちらは、anime_predict_dfのラベル列を使います。予測分布の推移のアニメーションを作成します。600)事後分布(ガンマ分布)の推移データが増えるにつれて、真のパラメータ付近の確率密度が大きくなっていくのを確認できます。ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.おわりに加筆修正の際に記事を分割しました。修正のペースを上げたい。2022.07.21:加筆修正しました。『パターン認識と機械学習』の独学時のまとめです。00:04何かおかしな内容があれば、ご指摘いただけるととても嬉しいです。【R】円周の作図5.3:正射影ベクトルの導出【『スタンフォード線形代数入門』のノート】スポンサーリンクはてなブログをはじめよう!anemptyarchiveさんは、はてなブログを使っています。ブログを報告する引用をストックしましたストック一覧を見る閉じる引用するにはまずログインしてくださいログイン閉じる引用をストックできませんでした。再度お試しください閉じる限定公開記事のため引用できません。