モンテカルロシミュレーションの実装-t統計量の真の有意水準-

著 J.アルバート「Rで学ぶベイズ統計学入門」のp11に基づいて作成

母集団が正規分布と等分散性の標準的な仮説に従わない場合、
t統計量の真の有意水準がどうなるか調べたいとしよう。
一般に真の有意水準は次の条項に依存する.

• あらかじめ定めた有意水準
• 母集団の形
• 二つの標準偏差ではかられる二つの母集団の広がり方
• 標本サイズmとn

次で表される真の有意水準を推定したい

を計算するシミュレーションアルゴリズムの概要を示す。

1. ランダム標本を第一の母集団から、を第２の母集団からシミュレーションする
2. 二つの標本からt統計量Tを計算する
3. |T|が棄却値を越えているかを判定し、越えていればを棄却する

アルゴリズムの三つのステップをN回繰り返す.
真の有意水準は以下で推定する.

= の棄却の数/N

以下はそのコード

#-*- coding: utf-8 -*-
from numpy import *
from scipy import stats

alpha=.1 #有意水準
m=10 #mの数
n=10 #nの数
N=1000.0 #シミュレーションの回数
reject=0 #棄却された回数を代入するオブジェクト

def tstatistic(x,y):
	m=len(x)
	n=len(y)
	#stats.nanstd(x)　で不偏分散を求める
	sp=sqrt(((m-1)*stats.nanstd(x)**2+(n-1)*stats.nanstd(y)**2)/(m+n-2))
	stat=(mean(x)-mean(y))/(sp*sqrt(1./m + 1./n))
	return stat

if __name__=="__main__":
	for i in range(1,N+1):
		x=random.normal(0,1,m)
		y=random.normal(0,1,n)
		a=tstatistic(x,y)
		if abs(a) > stats.t.ppf(1-alpha/2,n+m-2):
			reject=reject+1

	level=reject/N
	print  level