sage/graph/ギブス・サンプラー
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
単語検索
|
最終更新
|
ヘルプ
]
開始行:
[[FrontPage]]
#contents
2011/07/11からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://www15191ue.sakura.ne.jp:8000/home/pub/16/
また、Sageのサーバを公開しているサイト（http://www.sagenb...
アップロードし、実行したり、変更していろいろ動きを試すこ...
* Sageでグラフを再現してみよう：ギブス・サンプラー [#m9be...
この企画は、雑誌や教科書にでているグラフをSageで再現し、
グラフの意味を理解すると共にSageの使い方をマスターするこ...
今回は、
[[計算統計Ⅱマルコフ連鎖モンテカルロ法とその周辺>http://ww...
のp169のギブス・サンプラーの例（図２）を題材にします。
&ref(fig2.png);
** 図２について [#fb4ab664]
図２は、n=100個の乱数をN(5, 1)から発生させて、事前分布\(\...
\sigma \sim IG(0.001/2, 0.001/2)\)とおいてギブス・サンプ...
** 条件付き事後分布 [#wdeecc3b]
データを表現する条件付き事後分布からのサンプリングが可能...
ギブス・サンプリングを使うことができます。
平均\(\mu\)が分散\(\sigma^2\)の正規分布に従うとすると、事...
\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
$$
*** ギブス・サンプリングのアルゴリズム [#h22f6047]
事後分布からのサンプリングアルゴリズムは、以下のようにな...
+ 初期値(\(\mu^{(0)}, \sigma^{2(0)}\)を決め、t = 1とする
+ \(\mu\)の条件付き事後分布から、\(\mu^{(t)}\)をサンプリ...
$$
\mu^{(t)} | \sigma^{2(t-1)}, x \sim N(\mu_1^{(t)}, \sigma...
$$
$$
\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-...
{\sigma^{-2}_0 + n \sigma^{-2(t-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 -...
$$
+ t = t+1として、(2)に戻る
sageへの入力:
#pre{{
# Rのユーティリティ関数を読み込む
attach(DATA+'RUtil.py')
}}
*** 必要なライブラリの読み込み [#zd41c02d]
今回の例では、逆ガンマ分布のサンプリングが必要ですので、p...
使用することにしました。
sageへの入力:
#pre{{
# パッケージがインストールされていない場合に１度だけ実行
# r('install.packages("pscl")')
# igamma分布からサンプリングするためにライブラリpsclをロ...
r('library(pscl)')
}}
sageの出力:
#pre{{
[1] "pscl" "vcd" "colorspace" "grid" ...
[8] "lattice" "mvtnorm" "MASS" "stats" ...
[15] "datasets" "methods" "base"
}}
** テストデータの生成 [#dd51db26]
例題の説明から、正規分布N(5, 1)からサンプルを１００個生成...
通常なら、
#pre{{
X = [gauss(5,1) for i in range(100)]
}}
とするところですが、常に同じサンプルが生成できるようにRea...
(sage5.0から使えるようになった)を使ってテストデータを生成...
生成したデータをプロットしています。
sageへの入力:
#pre{{
# テストデータとして、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)
}}
&ref(out1.png);
*** テストデータの分布 [#le314dad]
生成したデータをヒスとグラムで表示すると以下のようになっ...
sageへの入力:
#pre{{
r(X).name('X')
graph = preGraph("hist.pdf")
r('hist(X)')
postGraph(graph)
}}
&ref(out2.png);
*** IG関数の定義 [#k90aa058]
\(\sigma^2\)の分布に逆ガンマ関数が使われているのは、その...
Sageには、逆ガンマ関数がないため、Rのrigamma関数を使って...
IG関数は以下のようになります。
sageへの入力:
#pre{{
# IGの定義
def IG(n, s):
return N(r('rigamma(1, %g, %g)'%(N(n), N(s))))
}}
** ギブス・サンプリングの初期設定 [#g0e5281c]
ギブス・サンプリングされた変数\(\mu, \sigma^2\)の値は、リ...
テストデータの個数n, xの平均をx_bar, \(\mu_0\)をmu_0、\(\...
\(n_0, S_0\)をそれぞれ、n_0, S_0変数にセットします。
\(sigma^2_0 = 1000\)としているのは、事前分布情報がない場...
表現しています。
また、\(n_0, S_0\)を0.001と小さく取ることによって１回のサ...
sageへの入力:
#pre{{
# 変数の初期設定
Sigma2= []
Mu = []
n = len(X)
x_bar = mean(X)
mu_0 = 0
sigma2_0 = 1000
n_0 = 0.001
S_0 = 0.001
}}
** ギブス・サンプリングの実行 [#qc428fec]
ギブス・サンプリングでは、最初の1000回までを稼働検査期間...
として、サンプリングから外しています。
その後の1000回をサンプリング期間とします。（モデルの複雑...
sageへの入力:
#pre{{
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 +...
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 +...
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))
}}
*** 結果のプロット [#y95d302d]
\(\mu, \sigma^2\)の変動と密度分布を以下にプロットします。
sageへの入力:
#pre{{
mu_plt = list_plot(Mu, plotjoined=True)
sig_plt = list_plot(Sigma2, plotjoined=True)
graphics_array([mu_plt, sig_plt]).show(figsize=[8, 3])
}}
&ref(out3.png);
sageへの入力:
#pre{{
# 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への入力:
#pre{{
html.table([[getGraph(mu_density, fac=0.5), getGraph(sig_...
}}
&ref(out4.png);
*** 標本の散布図 [#m12d88b7]
標本の散布図を以下に示します。サンプリングによって\(\mu, ...
sageへの入力:
#pre{{
list_plot(zip(Mu, Sigma2)).show(ymin=0.5, ymax=1.75, xmin...
}}
&ref(out5.png);
** コメント [#qf6eb026]
#vote(おもしろかった[3],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
- 逆ガンマ分布について、ベイズモデルでは正規分布の分散の...
#comment_kcaptcha
終了行:
[[FrontPage]]
#contents
2011/07/11からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://www15191ue.sakura.ne.jp:8000/home/pub/16/
また、Sageのサーバを公開しているサイト（http://www.sagenb...
アップロードし、実行したり、変更していろいろ動きを試すこ...
* Sageでグラフを再現してみよう：ギブス・サンプラー [#m9be...
この企画は、雑誌や教科書にでているグラフをSageで再現し、
グラフの意味を理解すると共にSageの使い方をマスターするこ...
今回は、
[[計算統計Ⅱマルコフ連鎖モンテカルロ法とその周辺>http://ww...
のp169のギブス・サンプラーの例（図２）を題材にします。
&ref(fig2.png);
** 図２について [#fb4ab664]
図２は、n=100個の乱数をN(5, 1)から発生させて、事前分布\(\...
\sigma \sim IG(0.001/2, 0.001/2)\)とおいてギブス・サンプ...
** 条件付き事後分布 [#wdeecc3b]
データを表現する条件付き事後分布からのサンプリングが可能...
ギブス・サンプリングを使うことができます。
平均\(\mu\)が分散\(\sigma^2\)の正規分布に従うとすると、事...
\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
$$
*** ギブス・サンプリングのアルゴリズム [#h22f6047]
事後分布からのサンプリングアルゴリズムは、以下のようにな...
+ 初期値(\(\mu^{(0)}, \sigma^{2(0)}\)を決め、t = 1とする
+ \(\mu\)の条件付き事後分布から、\(\mu^{(t)}\)をサンプリ...
$$
\mu^{(t)} | \sigma^{2(t-1)}, x \sim N(\mu_1^{(t)}, \sigma...
$$
$$
\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-...
{\sigma^{-2}_0 + n \sigma^{-2(t-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 -...
$$
+ t = t+1として、(2)に戻る
sageへの入力:
#pre{{
# Rのユーティリティ関数を読み込む
attach(DATA+'RUtil.py')
}}
*** 必要なライブラリの読み込み [#zd41c02d]
今回の例では、逆ガンマ分布のサンプリングが必要ですので、p...
使用することにしました。
sageへの入力:
#pre{{
# パッケージがインストールされていない場合に１度だけ実行
# r('install.packages("pscl")')
# igamma分布からサンプリングするためにライブラリpsclをロ...
r('library(pscl)')
}}
sageの出力:
#pre{{
[1] "pscl" "vcd" "colorspace" "grid" ...
[8] "lattice" "mvtnorm" "MASS" "stats" ...
[15] "datasets" "methods" "base"
}}
** テストデータの生成 [#dd51db26]
例題の説明から、正規分布N(5, 1)からサンプルを１００個生成...
通常なら、
#pre{{
X = [gauss(5,1) for i in range(100)]
}}
とするところですが、常に同じサンプルが生成できるようにRea...
(sage5.0から使えるようになった)を使ってテストデータを生成...
生成したデータをプロットしています。
sageへの入力:
#pre{{
# テストデータとして、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)
}}
&ref(out1.png);
*** テストデータの分布 [#le314dad]
生成したデータをヒスとグラムで表示すると以下のようになっ...
sageへの入力:
#pre{{
r(X).name('X')
graph = preGraph("hist.pdf")
r('hist(X)')
postGraph(graph)
}}
&ref(out2.png);
*** IG関数の定義 [#k90aa058]
\(\sigma^2\)の分布に逆ガンマ関数が使われているのは、その...
Sageには、逆ガンマ関数がないため、Rのrigamma関数を使って...
IG関数は以下のようになります。
sageへの入力:
#pre{{
# IGの定義
def IG(n, s):
return N(r('rigamma(1, %g, %g)'%(N(n), N(s))))
}}
** ギブス・サンプリングの初期設定 [#g0e5281c]
ギブス・サンプリングされた変数\(\mu, \sigma^2\)の値は、リ...
テストデータの個数n, xの平均をx_bar, \(\mu_0\)をmu_0、\(\...
\(n_0, S_0\)をそれぞれ、n_0, S_0変数にセットします。
\(sigma^2_0 = 1000\)としているのは、事前分布情報がない場...
表現しています。
また、\(n_0, S_0\)を0.001と小さく取ることによって１回のサ...
sageへの入力:
#pre{{
# 変数の初期設定
Sigma2= []
Mu = []
n = len(X)
x_bar = mean(X)
mu_0 = 0
sigma2_0 = 1000
n_0 = 0.001
S_0 = 0.001
}}
** ギブス・サンプリングの実行 [#qc428fec]
ギブス・サンプリングでは、最初の1000回までを稼働検査期間...
として、サンプリングから外しています。
その後の1000回をサンプリング期間とします。（モデルの複雑...
sageへの入力:
#pre{{
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 +...
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 +...
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))
}}
*** 結果のプロット [#y95d302d]
\(\mu, \sigma^2\)の変動と密度分布を以下にプロットします。
sageへの入力:
#pre{{
mu_plt = list_plot(Mu, plotjoined=True)
sig_plt = list_plot(Sigma2, plotjoined=True)
graphics_array([mu_plt, sig_plt]).show(figsize=[8, 3])
}}
&ref(out3.png);
sageへの入力:
#pre{{
# 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への入力:
#pre{{
html.table([[getGraph(mu_density, fac=0.5), getGraph(sig_...
}}
&ref(out4.png);
*** 標本の散布図 [#m12d88b7]
標本の散布図を以下に示します。サンプリングによって\(\mu, ...
sageへの入力:
#pre{{
list_plot(zip(Mu, Sigma2)).show(ymin=0.5, ymax=1.75, xmin...
}}
&ref(out5.png);
** コメント [#qf6eb026]
#vote(おもしろかった[3],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
- 逆ガンマ分布について、ベイズモデルでは正規分布の分散の...
#comment_kcaptcha
ページ名:
SmartDoc
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
