通常，データ分析では仮説検定を1回しか行わない．母集団から取り出したデータ（標本）が1つしかないからである．

しかしR（に類したプログラマブルな統計パッケージ）を利用すれば簡単に， 母集団からの標本抽出を人工的に繰り返すことができる． そこで仮説検定を繰り返して「p値の分布」を作り，1回限りの検定結果に基づいて 推論することの意味を確認してみよう.

最初に正規分布に従う乱数を作る.以下のコードは平均0，標準偏差1の正規分布に従う乱数を10個作る

rnorm(10, 0, 1)
##  [1] -0.458891671  0.806772309  0.851017773 -0.428373654  1.448784062
##  [6] -0.009394411 -2.536925047  0.755071620 -0.505671028 -1.124717739

これは直感的には確率密度関数が

curve(dnorm(x), -5, 5)

であるような正規分布を母集団として想定し，そこから10個の標本値を抽出したこと意味する． （あとで平均値の差の検定をするだけなので，母集団は正規分布でなくても構わない）

次の関数coinはn個の正規分布(平均0，標準偏差sd)に従う乱数を二列作り，その平均値に差があるかどうかをt検定する関数である

coin内の二つの分布の平均は同じだから，帰無仮説は正しい．

coin<-function(n,sd){
x <- rnorm(n, 0, sd)#n個の正規乱数をxに代入
y <- rnorm(n, 0, sd)#n個の正規乱数をyに代入
z<-t.test(x,y,var.equal=T)#
z$p.value}#t検定の結果からp値だけを返す
coin(1000,1)
## [1] 0.3553903

次の関数testは関数coinをk回反復する

test<-function(k,n,sd){
x <- c()   # 空の（要素数ゼロの）リストを作る
for (i in 1:k) {
  y <- coin(n,sd) #　n個のサンプルを使って，t検定の結果からp値のみ返す
  x <- c(x,  y)   }#結果ベクトルに追加していく
x}

p値のヒストグラムをプロットする

hist(test(k=1000,n=10000,sd=5), breaks = 100)

シミュレーションにより，帰無仮説が正しいときは，p値の分布が一様分布に近づく様子が分かる

次に関数を一般化して，帰無仮説が正しくない場合を表現しよう

coin2<-function(n,m1,m2,sd){
x <- rnorm(n, m1, sd)#n個の正規乱数をxに代入
y <- rnorm(n, m2, sd)#n個の正規乱数をyに代入
z<-t.test(x,y,var.equal=T)#
z$p.value}#t検定の結果からp値だけを返す

今度は帰無仮説が正しくない場合， m1,m2を明示的に指定して，母集団の平均値に差がある状態を作る．

coin2(n=100,m1=0,m2=1,sd=1)
## [1] 1.228741e-13
test2<-function(k,n,m1,m2,sd){
x <- c()   # 空の（要素数ゼロの）リストを作る
for (i in 1:k) {
  y <- coin2(n,m1,m2,sd) #　n個のサンプルを使って，t検定の結果からp値のみ返す
  x <- c(x,  y)   }#結果ベクトルに追加していく
x}

関数coin2を1000回繰り返して，t検定の p値の分布をプロットする

hist(test2(k=1000,n=10000,m1=0,m2=0.5,sd=5), breaks = 100)

対立仮説が正しい場合は，p値の分布が0付近に集中することが分かる 比較のため，平均値の差をさらに小さくしてみよう

hist(test2(k=1000,n=10000,m1=0,m2=0.1,sd=5), breaks = 100)

p値の分布が，変化する様子が分かる． 次に平均値を固定して，母集団の標準偏差だけをsd=5からsd=1変化させてみる．

hist(test2(k=1000,n=10000,m1=0,m2=0.1,sd=1), breaks = 100)

母分散が小さければ，小さなp値が出やすい．次に標本数だけ n=10000からn=100に変化させる

hist(test2(k=1000,n=100,m1=0,m2=0.1,sd=1), breaks = 100)

このように標本数が小さいと，帰無仮説を棄却できないケースが増える．