Tumblrはじめに『ベイズ推論による機械学習入門』の学習時のノートです。基本的な内容は「数式の行間を読んでみた」とそれを「RとPythonで組んでみた」になります。「数式」と「プログラム」から理解するのが目標です。この記事は、4.4.2項の内容です。)」、「事前分布をガウス・ウィシャート分布とディリクレ分布」とするガウス混合モデルに対するギブスサンプリングをPythonで実装します。省略してある内容等ありますので、本とあわせて読んでください。初学者な自分が理解できるレベルまで落として書き下していますので、分かる人にはかなりくどくなっています。参考文献おわりに・Pythonでやってみる人工的に生成したデータを用いて、ベイズ推論を行ってみましょう。アルゴリズム4.5を参考に実装します。利用するライブラリを読み込みます。pltscipyライブラリのstatsモジュールから、多次元ガウス分布multivariate_normal、ウィシャート分布wishart、ディリクレ分布dirichletのクラスを利用します。確率密度の計算にはpdf()、乱数の生成にはrvs()を使います。・観測モデルの設定まずは、観測モデルを設定します。\boldsymbol{\Lambda}_k^{-1})$とします。\boldsymbol{\pi})$によって行います。観測モデルのパラメータを設定します。2$のときのみ動作します。\sigma_{k,D,D}^2)$をsigma2_truth_ddkとして、値を指定します。$\boldsymbol{\Sigma}_k$は正定値行列です。\boldsymbol{\Sigma}_k^{-1}$を精度行列と呼びます。設定したパラメータは次のようになります。250.]]]混合比率を設定します。0$の値を指定します。1$となる)確率を表します。この3つのパラメータの値を観測データから求めるのがここでの目的です。設定した観測モデルをグラフで確認しましょう。作図用の点を作成します。x_{n,2})$がとり得る点(値の組み合わせ)を作成します。$x_{n,1}$がとり得る値(x軸の値)をnp.linspace()で作成してx_1_lineとします。この例では、$K$個のクラスタそれぞれで平均値から標準偏差の3倍を引いた値と足した値を計算して、その最小値から最大値までを範囲とします。np.linspace()を使うと指定した要素数で等間隔に切り分けます。np.arange()を使うと切り分ける間隔を指定できます。処理が重い場合は、この値を調整してください。同様に、$x_{n,2}$がとり得る値(y軸の値)も作成してx_2_lineとします。x_2_gridとします。これは、確率密度を等高線図にする際に格子状の点(2軸の全ての値が直交する点)を渡す必要があるためです。x_2_gridの要素を列として持つ2次元配列を作成してx_pointとします。これは、確率密度の計算に使います。計算結果の確率密度は1次元配列になります。x_2_gridと同じ形状にする必要があるため、形状をx_dimとして保存しておきます。作成した$\mathbf{x}$の点は次のようになります。]]1列目がx軸の値、2列目がy軸の値に対応しています。観測モデルの確率密度を計算します。sigmga2_truth_ddkを使って、x_point_matの点ごとに、次の混合ガウス分布の式で確率密度を計算します。$$$K$個の多次元ガウス分布に対して、$\boldsymbol{\pi}$を使って加重平均をとります。多次元ガウス分布の確率密度は、SciPyライブラリのmultivariate_normal.pdf()で計算できます。$k$番目の分布の場合は、データの引数Xにx_point、平均の引数muにmu_truth_kd[k]、分散共分散行列の引数sigmaにsigma2_truth_kdd[k]を指定します。この計算をfor文でK回繰り返し行います。各パラメータによって求めた$K$回分の確率密度をそれぞれpi_truth_k[k]で割り引いてmodel_probに加算していきます。計算結果は次のようになります。1.77630256e-12]観測モデルを作図します。等高線グラフはplt.contour()で描画します。$K$個の分布を重ねたような分布になります。積分して1となるようにしているため、$K$個の分布そのものよりも1つ1つの山が小さく(確率密度が小さく)なっています。続いて、構築した観測モデルに従ってデータを生成します。s_{n,K})$を生成し、各データにクラスタを割り当てます。x_{n,D})$を生成します。$N$個の潜在変数$\mathbf{S}$を生成します。size=N)生成するデータ数$N$をNとして値を指定します。カテゴリ分布に従う乱数は、多項分布に従う乱数生成関数np.random.multinomial()のn引数を1にすることで生成できます。確率(パラメータ)の引数probにpi_truth_k、試行回数の引数sizeにNを指定します。生成した$N$個の$K$次元ベクトルをs_truth_nkとします。結果は次のようになります。1]]これが本来は観測できない潜在変数であり、真のクラスタです。各行がデータを表し、値が1の列番号が割り当てられたクラスタ番号を表します。s_truth_nkから各データのクラスタ番号を抽出します。2]np.where()を使って行(データ)ごとに要素が1の列番号を抽出します。各データに割り当てられたクラスタに従い$N$個のデータを生成します。\boldsymbol{\Sigma}_k$を持つ多次元ガウス分布に従いデータ(値)を生成します。多次元ガウス分布の乱数は、np.random.multivariate_normaln()で生成できます。データごとにクラスタが異なるので、for文(リスト内包表記)で1データずつ処理します。$n$番目のデータのクラスタ番号s_true_n[n]を取り出してkとします。そのクラスタに従ってデータ$\mathbf{x}_n$を生成します。平均の引数meanにmu_true_kd[k]、分散共分散行列の引数covにsigma2_true_kdd[k]を指定します。1データずつ生成するので、データ数の引数sizeは1です。2次元配列で出力されるため、flatten()で1次元配列に変換しておきます。生成した観測データは次のようになります。-1.21215116e+01]]観測データの散布図を観測モデルに重ねて作図します。散布図はplt.scatter()で描画します。k)で抽出してk_idxとします。k_idxを使ってクラスタ$k$が割り当てられたデータを取り出して、クラスタごとに色を変えた散布図を描きます。観測データから各データのクラスタも推定します。\boldsymbol{\pi})$に対する共役事前分布を設定します。\boldsymbol{\alpha})$を設定します。各事前分布のパラメータ(超パラメータ)を設定します。m_D)$をm_d、精度行列パラメータの係数$\beta$をbetaとして値を指定します。w_{D,D})$をw_dd、自由度$\nu$をnuとして値を指定します。0$の値を指定します。パラメータを設定できたので、$\boldsymbol{\mu}$の事前分布をグラフで確認しましょう。作図用の$\boldsymbol{\mu}$の点を作成します。300)$\mathbf{x}$の点のときと同様に作成します。各軸の描画範囲は、真の平均の最小値から事前分布の標準偏差を引いた値から、真の平均の最大値に標準偏差を足した値までとします。\boldsymbol{\Lambda}_k)^{-1}$の対角成分の平方根です。この時点では$\boldsymbol{\Lambda}$が未知の値なので、代わりに$\boldsymbol{\Lambda}$の期待値$\mathbb{E}[\boldsymbol{\Lambda}]$を使うことにします。逆行列の計算はnp.linalg.inv()、対角成分の抽出はdiagonal()、平方根の計算はnp.sqrt()で行えます。面倒でしたら適当な値を直接指定してください。$\boldsymbol{\mu}$の事前分布の確率密度を計算します。w_dd))観測モデルのときと同様に、多次元ガウス分布の確率密度を計算します。\mathbb{E}[\boldsymbol{\Lambda}])^{-1}$とします。計算結果は次のようになります。1.19655008e-05]$\boldsymbol{\mu}$の事前分布を作図します。真の平均値をプロットしておきます。この真の平均を分布推定します。$\boldsymbol{\Lambda}$と$\boldsymbol{\pi}$の事前分布(と事後分布)については省略します。項「カテゴリ分布の学習と予測」と3.4.2項「多次元ガウス分布の学習と予測:精度が未知の場合」を参照してください。\boldsymbol{\pi}$を生成して初期値とします。\boldsymbol{\Lambda}$を生成します。).flatten()精度行列パラメータ$\boldsymbol{\Lambda}$の事前分布(ウィシャート分布)と平均パラメータ$\boldsymbol{\mu}$の事前分布(多次元ガウス分布)からそれぞれ$K$個のパラメータをサンプルして初期値とします。$$ウィシャート分布に従う乱数は、SicPyライブラリのwishart.rvs()で生成できます。クラスタごとに生成するので、データ数の引数sizeに1を指定します。自由度の引数dfにnu、スケール行列の引数scaleにw_ddを指定します。多次元ガウス分布の乱数生成関数np.random.multivariate_normal()の引数meanにm_d、covには今サンプルしたlambda_kdd[k]と係数betaの積をnp.linalg.inv()で逆行列に変換して指定します。2次元配列で出力されるため、flatten()で1次元配列に変換しておきます。生成したパラメータを確認します。2.22053690e-03]]]$\boldsymbol{\pi}$の事前分布から$\boldsymbol{\pi}$を生成します。size=1).flatten()$\boldsymbol{\pi}$の事前分布(ディリクレ分布)に従い$\boldsymbol{\pi}$をサンプルして初期値とします。dirichlet.rvs()で生成できます。パラメータの引数alphaにalpha_kを指定します。2次元配列で出力されるため、flatten()で1次元配列に変換しておきます。生成したパラメータを確認します。0.82242016]1.0$K$個の要素の和が1になるのを確認できます。パラメータの初期値が得られたので、グラフで確認しましょう。初期値を使った分布を計算します。観測モデルのときと同様に作図します。サンプルを繰り返すことで真の分布に近付けていきます。以上でモデルが整いました。次からは、推論を行っていきます。\boldsymbol{\pi}$のサンプリングを交互に行います。\boldsymbol{\pi})$を計算します。得られた条件付き分布から$\mathbf{S}$をサンプルします。\mathbf{S})$を計算します。得られた条件付き分布からパラメータをサンプルします。ここまでが1回の試行です。次の試行では、前回サンプルしたパラメータを用いて$\mathbf{S}$の条件付き分布を更新となります。この処理を指定した回数繰り返します。各パラメータの変数に関して、添字を使って代入するため予め変数を用意しておく必要があります。後で各パラメータの更新値や分布の推移を確認するために、それぞれ試行ごとにリスト型のtrace_***に値(NumPy配列)を保存していきます。関心のあるパラメータを記録してください。'%)')・処理の解説ギブスサンプリングで行う処理を確認していきます。・潜在変数のサンプリング潜在変数$\mathbf{S}$の条件付き分布(のパラメータ)を計算して、新たな$\mathbf{S}$をサンプルします。$\mathbf{S}$の条件付き分布のパラメータを計算します。\eta_{n,K})$の各項は、次の式で計算します。1$となる)確率を表します。\boldsymbol{\mu}_k)$の計算を1から$N$の全ての組み合わせで行います(tmp_eta_nの計算)。N$のマトリクスになります。\boldsymbol{\mu}_k)$の計算結果です。これは不要なので、np.diag()で対角成分(同じデータによる計算結果)のみを取り出します。0.0000001$)を加えています。1.]行(データ)ごとの和が1になるのを確認できます。eta_nkを用いて、新たな$\mathbf{S}$を生成します。size=1).flatten()真のクラスタのときと同様に、潜在変数をサンプルします。$$ただし、ここではデータによってパラメータ(各クラスタの割り当て確率)が異なるため、for文で1データずつ生成します。0.]]これで$\mathbf{S}$が得られました。\boldsymbol{\pi}$を更新します。\boldsymbol{\Lambda}$をサンプルします。$\boldsymbol{\mu}_k$の条件付き分布のパラメータを計算します。$$で計算して、それぞれbeta_hat_k、m_hat_kdに代入します。-33.70478007]対象としているクラスタのデータには1を掛けてそのままの値で、それ以外のクラスタのデータには0を掛けて値が0になっています。対象となるクラスタのデータのみの和をとっています。s_{n,k}$は各クラスタが割り当てられたデータ数を表します。$\boldsymbol{\lambda}_k$の条件付き分布のパラメータを計算します。$$で計算して、結果をそれぞれnu_hat_k、w_hat_ddkに代入します。w_hat_ddkを用いて、新たな$\boldsymbol{\lambda}_k$を生成します。scale=w_hat_kdd[k])パラメータ$\hat{\mathbf{W}}_k$を持つ自由度$\hat{\nu}_k$のウィシャート分布から$\boldsymbol{\lambda}_k$をサンプルします。$$初期値のときと同様に、wishart.rvs()の自由度の引数dfにnu_hat_k[k]、スケール行列の引数scaleにw_hat_kdd[k]を指定します。beta_hat_kとlambda_ddkを用いて、新たな$\boldsymbol{\mu}_k$を生成します。\boldsymbol{\lambda}_k$の多次元ガウス分布から$\boldsymbol{\mu}_k$をサンプルします。$$np.random.multivariate_normal()の平均の引数meanにm_hat_kd[k]、分散共分散行列の引数covに今サンプルしたlambda_kdd[k]と係数beta_hat_k[k]との積をnp.linalg.inv()で逆行列に変換して指定します。\hat{\mathbf{W}}_K\}$が得られます。・混合比率のサンプリング混合比率パラメータ$\boldsymbol{\pi}$の条件付き分布(のパラメータ)を計算し、新たな$\boldsymbol{\pi}$をサンプルします。$\boldsymbol{\pi}$の条件付き分布のパラメータを計算します。$$で計算して、結果をalpha_hat_kとします。120.]alpha_hat_kを用いて、新たな$\boldsymbol{\pi}$を生成します。alpha=alpha_hat_k).flatten()パラメータ$\hat{\boldsymbol{\alpha}}$を持つディリクレ分布から$\boldsymbol{\pi}$をサンプルします。$$事前分布からのサンプルのときと同様に、dirichlet.rvs()の引数alphaにalpha_hat_kを指定します。0.47247398]1.0以上がギブスサンプリングによる近似推論で行うの個々の処理です。これを指定した回数繰り返し行うことで、各パラメータの値を徐々に真の値に近付けていきます。推論結果を確認していきます。beta_hat_kとlambda_ddkの最後の更新値を用いて、$\boldsymbol{\mu}$の事後分布を計算します。)事前分布のときと同様に計算します。超パラメータが$K$個になったので、多次元ガウス分布も$K$個になります。$K$回分の計算結果を2次元配列に格納していきます。\boldsymbol{\Lambda}_k$です。作成した配列は次のようになります。0.00000000e+000]]$\boldsymbol{\mu}$の事後分布を作図します。事前分布のときと同様に作図できます。$K$個のガウス分布がそれぞれ真の平均付近をピークとする分布を推定できています。lambda_ddkとpi_kを用いて、混合分布を計算します。上はplt.contour()、下はplt.contourf()で作成したものです。クラスタリングできていることが確認できます。クラスタの番号についてはランダムに決まるため真のクラスタの番号と一致しません。\hat{\boldsymbol{\alpha}}$の更新値の推移を確認します。・コード(クリックで展開)学習ごとの値を格納した各パラメータのリストtrace_***をnp.arange()でNumPy配列に変換して作図します。\boldsymbol{\pi}$の値の推移を確認しましょう。・コード(クリックで展開)同様に作図します。それぞれ真の値***_truthを水平線plt.hlines()で描画しておきます。試行回数が増えるに従って真の値に近づいているのを確認できます。・おまけ:アニメーションによる推移の確認animationモジュールを利用して、各分布の推移のアニメーション(gif画像)を作成するためのコードです。上はplt.contour()、下はplt.contourf()で作成したものです。(アップロードサイズの上限を超えたためアニメーションを分割しました。左が0から100回目の結果まで、右が101から200回目の結果です。ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.おわりに今回は先にPythonで組んでみました。Rではどこかミスっててまだうまく推定できてません。この記事の投稿日の朝に公開されたモーニング娘。LOVE」です!歯磨けよ~【次節の内容】www.anarchive-beta.com2021/05/14:加筆修正しました。『パターン認識と機械学習』の独学時のまとめです。00:04何かおかしな内容があれば、ご指摘いただけるととても嬉しいです。【R】円周の作図5.3:正射影ベクトルの導出【『スタンフォード線形代数入門』のノート】スポンサーリンクはてなブログをはじめよう!anemptyarchiveさんは、はてなブログを使っています。ブログを報告する引用をストックしましたストック一覧を見る閉じる引用するにはまずログインしてくださいログイン閉じる引用をストックできませんでした。再度お試しください閉じる限定公開記事のため引用できません。