graph_gibbs_sampling

1383 days ago by takepwave

Hiroshi TAKEMOTO (take@pwv.co.jp)

Sageでグラフを再現してみよう:ギブス・サンプラー

この企画は、雑誌や教科書にでているグラフをSageで再現し、 グラフの意味を理解すると共にSageの使い方をマスターすることを目的としています。

今回は、計算統計Ⅱマルコフ連鎖モンテカルロ法とその周辺 のp169のギブス・サンプラーの例(図2)を題材にします。

図2について

図2は、n=100個の乱数をN(5, 1)から発生させて、事前分布$\mu \sim N(0, 1000), \sigma \sim IG(0.001/2, 0.001/2)$とおいてギブス・サンプリングを行った結果です。

条件付き事後分布

データが既存分布(今回の場合正規分布)を示している場合、ギブス・サンプリングを使って、 条件付き事後分布からサンプリングが可能になります。

平均$\mu$が分散$\sigma^2$の正規分布に従うとすると、事前分布$ \mu \sim N(\mu_0, \sigma^2_0), \sigma^2 \sim IG(n_0/2, S_0/2)$とすると、 $\mu$の条件付き事後分布は、以下のようになります。 $$ \mu | \sigma^2, x \sim N(\mu_1, \sigma^2_1), $$ $$ \sigma^{-2}_1 = \sigma^{-2}_0 + n \sigma^{-2}, \mu_1 = \frac{\sigma^{-2}_0 \mu_0 + n \sigma^{-2}\bar{x}} {\sigma^{-2}_0 + n \sigma^{-2}} $$

また、$\sigma^2$の条件付き事後分布は、以下のようになります。 $$ \sigma^2 | \mu, x \sim IG(n_1/2, S_1/2), $$ $$ n_1 = n_0 + n, S_1 = S_0 + \sum^{n}_{i = 1} ( x_i - \mu)^2 $$

ギブス・サンプリングのアルゴリズム

事後分布からのサンプリングアルゴリズムは、以下のようになります。

  1. 初期値($\mu^{(0)}, \sigma^{2(0)}$を決め、t = 1とする
  2. $\mu$の条件付き事後分布から、$\mu^{(t)}$をサンプリングする
  3. $$ \mu^{(t)} | \sigma^{2(t-1)}, x \sim N(\mu_1^{(t)}, \sigma^{2(t)}_1), $$ $$ \sigma^{-2(t)}_1 = \sigma^{-2}_0 + n \sigma^{-2(t-1)}, \mu_1^{(t)} = \frac{\sigma^{-2}_0 \mu_0 + n \sigma^{-2(t-1)}\bar{x}} {\sigma^{-2}_0 + n \sigma^{-2(t-1)}} $$
  4. $\sigma$の条件付き事後分布から、$\sigma^{2(t)}$をサンプリングする
  5. $$ \sigma^{2(t)} | \mu^{(t)}, x \sim IG(n_1/2, S_1^{(t)}/2), $$ $$ n_1 = n_0 + n, S_1^{(t)} = S_0 + \sum^{n}_{i = 1} ( x_i - \mu^{(t)})^2 $$
  6. t = t+1として、(2)に戻る

# RとPandasのデータフレームを相互に変換する関数を読み込む # Rの必要なライブラリ r('library(ggplot2)') r('library(jsonlite)') # RUtilの読み込み load(DATA + 'RUtil.py') 
       

必要なライブラリの読み込み

今回の例では、逆ガンマ分布のサンプリングが必要ですので、psclパッケージのigamma関数を 使用することにしました。

# パッケージがインストールされていない場合に1度だけ実行 # r('install.packages("pscl")') # igamma分布からサンプリングするためにライブラリpsclをロード r('library(pscl)') 
       
 [1] "pscl"      "lattice"   "MASS"      "jsonlite"  "ggplot2"   "stats"
"graphics"  "grDevices"
 [9] "utils"     "datasets"  "methods"   "base"     
 [1] "pscl"      "lattice"   "MASS"      "jsonlite"  "ggplot2"   "stats"     "graphics"  "grDevices"
 [9] "utils"     "datasets"  "methods"   "base"     

テストデータの生成

例題の説明から、正規分布N(5, 1)からサンプルを100個生成します。 通常なら、

	X = [gauss(5,1) for i in range(100)]
とするところですが、常に同じサンプルが生成できるようにRealDistribution関数 (sage5.0から使えるようになった)を使ってテストデータを生成します。

生成したデータをプロットしています。

# テストデータとして、N(5,1) #X = [gauss(5,1) for i in range(100)] # 必ず同じ分布になるようにRealDistributionに変更 T = RealDistribution('gaussian', 1, seed=100) X = [T.get_random_element()+5 for i in range(100)] # サンプルのプロット list_plot(X).show(figsize=5) 
       

テストデータの分布

生成したデータをヒスとグラムで表示すると以下のようになっています。

r(X).name('X') graph = preGraph("hist.pdf") r('hist(X)') postGraph(graph) 
       

IG関数の定義

$\sigma^2$の分布に逆ガンマ関数が使われているのは、その共役分布が正規分布であることに起因しています。

Sageには、逆ガンマ関数がないため、Rのrigamma関数を使って逆ガンマ分布のサンプリングをします。 IG関数は以下のようになります。

# IGの定義 def IG(n, s): return sageobj(r('rigamma(1, %g, %g)'%(N(n), N(s)))) 
       

ギブス・サンプリングの初期設定

ギブス・サンプリングされた変数$\mu, \sigma^2$の値は、リストMu, Sigma2にセットします。

テストデータの個数n, xの平均をx_bar, $\mu_0$をmu_0、$\sigma^2_0$をsigma2_0に $n_0, S_0$をそれぞれ、n_0, S_0変数にセットします。

$sigma^2_0 = 1000$としているのは、事前分布情報がない場合、分散を大きく取って不確実性を 表現しています。

また、$n_0, S_0$を0.001と小さく取ることによって1回のサンプリングによる変動を小さくしています。

# 変数の初期設定 Sigma2= [] Mu = [] n = len(X) x_bar = mean(X) mu_0 = 0 sigma2_0 = 1000 n_0 = 0.001 S_0 = 0.001 
       

ギブス・サンプリングの実行

ギブス・サンプリングでは、最初の1000回までを稼働検査期間(buring in period) として、サンプリングから外しています。

その後の1000回をサンプリング期間とします。(モデルの複雑さなどで回数は異なります)

sigma2 = sigma2_0 # 稼働検査期間(buring in period) for i in range(1000): # muの更新 sigma2_1 = 1/(1/sigma2_0 + n /sigma2) mu_1 = (mu_0/sigma2_0 + n/sigma2*x_bar)/(1/sigma2_0 + n/sigma2) mu = gauss(mu_1, sqrt(sigma2_1)) # sigmaの更新 n_1 = n_0 + n S_1 = S_0 + sum([(x - mu)^2 for x in X]) sigma2 = IG(n_1/2, S_1/2) # サンプリング for i in range(1000): # muの更新 sigma2_1 = 1/(1/sigma2_0 + n /sigma2) mu_1 = (mu_0/sigma2_0 + n/sigma2*x_bar)/(1/sigma2_0 + n/sigma2) mu = gauss(mu_1, sqrt(sigma2_1)) #sigmaの更新 n_1 = n_0 + n S_1 = S_0 + sum([(x - mu)^2 for x in X]) sigma2 = IG(n_1/2, S_1/2) # サンプリングを追加 Mu.append(mu) Sigma2.append(N(sigma2)) 
       

結果のプロット

$\mu, \sigma^2$の変動と密度分布を以下にプロットします。

mu_plt = list_plot(Mu, plotjoined=True) sig_plt = list_plot(Sigma2, plotjoined=True) graphics_array([mu_plt, sig_plt]).show(figsize=[8, 3]) 
       
# Muの密度分布図の作成 r(Mu).name('Mu') mu_density = preGraph("mu_density.pdf") r('plot(density(Mu))') offGraph() # Sigma2の密度分布図の作成 r(Sigma2).name('Sigma2') sig_density = preGraph("sig_density.pdf") r('plot(density(Sigma2))') offGraph() 
       
html.table([[getGraph(mu_density, fac=0.5), getGraph(sig_density, fac=0.5)]]) 
       

標本の散布図

標本の散布図を以下に示します。サンプリングによって$\mu, \sigma^2$が不変分布に収束していることが分かります。

list_plot(zip(Mu, Sigma2)).show(ymin=0.5, ymax=1.75, xmin=4.4, xmax=5.5, figsize=5)