1. No category

数値シミュレーション 第 15 回 教材
担 当 神保 雅人
2014 年 1 月 24 日（金）実施
ベイズ統計 3
自然共役事前分布の利用 (続き)
今回も，母集団が正規分布の場合のベイズ推測を取り扱う。
【平均が既知の場合】
母集団の平均μの値(μ=c )は既知，分散σ2 の値は未知とする。このとき，n 個のデータ
x
( x1 , x2 ,
, xn ) を得た場合，このデータの母集団の平均に対する分散 s 2
1
n
(標本分散)に基づいて，尤度関数は
2
f (x |
n
)
e
1
2
i 1
1
2 2
( xi c ) 2
2
(
n
)
2
n
( xi c) 2
i 1
ns 2
e
2 2
となる。この場合，σ2 の事前分布を逆ガンマ分布 IG (n0/2, n0S0/2)，即ち
(
2
)
(
2
)
n 0 S0
n0
( 2 1)
e
2 2
とすると，事後分布も逆ガンマ分布 IG ([n+n0]/2, [ns 2+n0S0]/2) となる。
【平均及び分散が共に未知の場合】
平均μ，分散σ2 の値が共に未知とする。このとき，n 個のデータ x
場合，尤度関数は，x
2
f (x | ,
)
1
n
n
xi ，
n 1 ，s
2
1
i 1
n
( x1 , x2 ,
, xn ) を得た
x ) 2 として
( xi
i 1
e
1
2
i 1
n
1
2 2
( xi
)2
(
2
)
n
2
e
1
2 2
[ s 2 n(
x )2 ]
となる。この場合，μ及びσ2 の同時自然共役事前分布は
( |
2
) (
2
)
(
2
)
n 1
( 02 1)
e
1 [n s
0 0 m0 (
2 2
0)
2
]
の様に，N(μ0,σ2/m0 )の分布と IG (n0/2, n0S0/2) の分布との積となる。
例題 1
自然共役事前分布から事後分布を求める次のプログラムを入力，ビルドして，実行せよ。ここ
で，クラス名は Sample15_1，ソースファイル名は Sample15_1.java とする。また，テキスト形式
－ 1 －
数値シミュレーション 第 15 回 教材
担 当 神保 雅人
のファイルを GNUPLOT で読み込み，plot コマンドにより，グラフを描画せよ。その際，ファ
イル名の指定の後に with linespoints のオプションを加える。
import
import
import
import
import
java.io.BufferedWriter;
java.io.FileNotFoundException;
java.io.FileWriter;
java.io.IOException;
java.io.PrintWriter;
public class Sample15_1 {
public static void main(String[] args) {
final int M=100,n=10,n0=10;
final double s2=0.1,S0=1.5,MAX=2.2;
double sig2=0.2,div=(MAX-sig2)/M;
String fname = "bayesian2.txt";
try {
PrintWriter pw = new PrintWriter(new BufferedWriter(new FileWriter(fname)));
for (int i=0;i<=M;i++) {
pw.printf("%f %f%n", sig2, IG(sig2,n/2-1,n*s2/2)*IG(sig2,n0/2,n0*S0/2));
sig2 += div;
}
pw.close();
} catch (FileNotFoundException fe) {
System.out.println(fname +"というファイルが開けません");
} catch (IOException err) {
System.out.println("IOException をキャッチ"+err);
}
}
private static double IG(double sig2,double a,double b) {
return Math.pow(sig2,-(a+1))*Math.exp(-b/sig2);
}
}
＊ このプログラムでは，尤度関数と事前分布とを掛け合わせて事後分布を求めているが，それぞ
れの正規化の因子を掛けていない。
例題 2
自然共役事前分布から事後分布を求める次のプログラムを入力，ビルドして，実行せよ。ここ
で，クラス名は Sample15_2，ソースファイル名は Sample15_2.java とする。また，テキスト形式
のファイルを GNUPLOT で読み込み，splot コマンドにより，グラフを描画せよ。その際，ファ
－ 2 －
数値シミュレーション 第 15 回 教材
担 当 神保 雅人
イル名の指定の後に with linespoints のオプションを加える。
import
import
import
import
import
java.io.BufferedWriter;
java.io.FileNotFoundException;
java.io.FileWriter;
java.io.IOException;
java.io.PrintWriter;
public class Sample15_2 {
public static void main(String[] args) {
final int M=20,n=10,nu=n-1,n0=10,m0=2;
final double s2=0.1,xave=1.0,S0=1.5,mu0=1.2,MAXm=1.5,MAXs=2.2,sig2init=0.2;
double mu=0.1,sig2=sig2init,divm=(MAXm-mu)/M,divs=(MAXs-sig2)/M;
String fname = "bayesian3.txt";
try {
PrintWriter pw = new PrintWriter(new BufferedWriter(new FileWriter(fname)));
double f,Pb;
for (int i=0;i<=M;i++) {
for (int j=0; j<=M; j++) {
f = N(mu,xave,sig2/n)*IG(sig2,n/2-1,nu*s2/2);
Pb = N(mu,mu0,sig2/m0)*IG(sig2,n0/2,n0*S0/2);
pw.printf("%f %f %E%n", sig2,mu,f*Pb);
sig2 += divs;
}
pw.println();
sig2=sig2init;
mu += divm;
}
pw.close();
} catch (FileNotFoundException fe) {
System.out.println(fname +"というファイルが開けません");
} catch (IOException err) {
System.out.println("IOException をキャッチ"+err);
}
}
private static double N(double mu,double ave,double sig2) {
return Math.exp(-(ave-mu)*(ave-mu)/(2.0*sig2));
}
private static double IG(double sig2,double a,double b) {
return Math.pow(sig2,-(a+1))*Math.exp(-b/sig2);
}
－ 3 －
数値シミュレーション 第 15 回 教材
担 当 神保 雅人
}
＊ このプログラムでは，尤度関数と事前分布とを掛け合わせて事後分布を求めているが，それぞ
れの正規化の因子を掛けていない。
演習 1
例題 1 のプログラムを参考にして，σ2 の値 σ2 に対する事前分布の値を，MAX の値は 5.2 に変
更して，テキスト形式のファイル bayesian2ex.txt に書き出す様に書き換えた Java プログラムを
作成し，実行せよ。ここで，クラス名は Ex15_1，ソースファイル名は Ex15_1.java とする。また，
ここで得られたテキスト形式のファイル及び例題で得られたテキスト形式のファイルを
GNUPLOT で読み込み，plot コマンドにより，両者を重ねたグラフを描画せよ。その際，ファイ
ル名の指定の後に with linespoints のオプションを加える。
演習 2
例題 2 のプログラムを参考にして，σ2 の値 μの値 σ2 ，μに対する事前分布の値を，MAXm
の値は 3.0 に変更して，
テキスト形式のファイル bayesian3ex.txt に書き出す様に書き換えた Java
プログラムを作成し，実行せよ。ここで，クラス名は Ex15_2，ソースファイル名は Ex15_2.java と
する。また，ここで得られたテキスト形式のファイル及び例題で得られたテキスト形式のファイ
ルを GNUPLOT で読み込み，splot コマンドにより，両者を重ねたグラフを描画せよ。その際，
ファイル名の指定の後に with linespoints のオプションを加える。
今回の提出物
1) 演習 1 で GNUPLOT を用いて作成したグラフのウィンドウで
『Toggle grid』
ボタンを押して，
グリッドを表示させたものを『Copy the plot to clipboard』ボタンを押してクリップボードにコ
ピーし，ペイントに貼り付けて，Ex15_1.jpg の名前で保存したもの(保存時のファイルの種類に注
意)
2) 演習 2 で GNUPLOT を用いて作成したグラフのウィンドウで
『Toggle grid』
ボタンを押して，
グリッドを表示させたものを『Copy the plot to clipboard』ボタンを押してクリップボードにコ
ピーし，ペイントに貼り付けて，Ex15_2.jpg の名前で保存したもの(保存時のファイルの種類に注
意)
－ 4 －