FrontPage

2011/07/11からのアクセス回数 6346

ここで紹介したSageワークシートは、以下のURLからダウンロードできます。

http://www15191ue.sakura.ne.jp:8000/home/pub/16/

また、Sageのサーバを公開しているサイト（http://www.sagenb.org/, http://www15191ue.sakura.ne.jp:8000/）にユーザIDを作成することで、ダウンロードしたワークシートを アップロードし、実行したり、変更していろいろ動きを試すことができます。

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

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

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

fig2.png

↑

図２について †

図２は、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とする
  1. \(\mu\)の条件付き事後分布から、\(\mu^{(t)}\)をサンプリングする $$ \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)}} $$
  1. \(\sigma\)の条件付き事後分布から、\(\sigma^{2(t)}\)をサンプリングする $$ \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 $$
  1. t = t+1として、(2)に戻る

sageへの入力:

# Rのユーティリティ関数を読み込む
attach(DATA+'RUtil.py')
↑

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

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

sageへの入力:

# パッケージがインストールされていない場合に１度だけ実行
# r('install.packages("pscl")')
# igamma分布からサンプリングするためにライブラリpsclをロード
r('library(pscl)')

sageの出力:

 [1] "pscl"       "vcd"        "colorspace" "grid"       "gam"        "splines"    "coda"      
 [8] "lattice"    "mvtnorm"    "MASS"       "stats"      "graphics"   "grDevices"  "utils"     
[15] "datasets"   "methods"    "base"
↑

テストデータの生成 †

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

	X = [gauss(5,1) for i in range(100)]

とするところですが、常に同じサンプルが生成できるようにRealDistribution関数 (sage5.0から使えるようになった)を使ってテストデータを生成します。

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

sageへの入力:

# テストデータとして、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)

out1.png

↑

テストデータの分布 †

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

sageへの入力:

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

out2.png

↑

IG関数の定義 †

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

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

sageへの入力:

# IGの定義
def IG(n, s):
    return N(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と小さく取ることによって１回のサンプリングによる変動を小さくしています。

sageへの入力:

# 変数の初期設定
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回をサンプリング期間とします。（モデルの複雑さなどで回数は異なります）

sageへの入力:

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\)の変動と密度分布を以下にプロットします。

sageへの入力:

mu_plt = list_plot(Mu, plotjoined=True)
sig_plt = list_plot(Sigma2, plotjoined=True)
graphics_array([mu_plt, sig_plt]).show(figsize=[8, 3])

out3.png

sageへの入力:

# 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()

sageへの入力:

html.table([[getGraph(mu_density, fac=0.5), getGraph(sig_density, fac=0.5)]])

out4.png

↑

標本の散布図 †

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

sageへの入力:

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

out5.png

↑

コメント †

選択肢 投票
おもしろかった 1  
そうでもない 0  
わかりずらい 0  

皆様のご意見、ご希望をお待ちしております。


(Input image string)
トップ   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
SmartDoc