カイ二乗分布とt分布

1 カイ二乗分布

標準正規分布の平方和はカイ二乗分布に従うことが知られています。

1.1 母平均が既知のとき

\(X\sim N(0,1)\) は標準正規分布に従う確率変数で、Xから得られたデータ \((x_1,x_2,\cdots,x_n)\) は独立に標準正規分布に従うとします。
独立な標準正規分布の偏差平方和はカイ二乗値であり、
　　　\(\chi^2=x_1^2+x_2^2+\cdots+x_n^2\)

標本分散は、
　　　\(V=\dfrac{x_1^2+x_2^2+\cdots+x_n^2}{n}=\dfrac{\chi^2}{n}\)

分散は カイ二乗値をnで割ったものです。（偏差平方の平均が分散ということになります）

n=10 の場合でシミュレーションしてみます。 （10個の和は自由度10）

X = c()
for (i in 1:10000) {
  x <- rnorm(10,0,1)
  s <- sum((x-0)^2)
  X <-c(X,s)
}
mean(X)

## [1] 9.959046

var(X)

## [1] 19.68336

hist(X,prob=T)
curve(dchisq(x,10),add=T)

n=10 の場合、 　　期待値 \(E(\chi^2)=10\) 、分散 \(V(\chi^2)=20\)
理論上の確率密度関数とも合致しているようです。

n個の平方和 \(\chi^2=x_1^2+x_2^2+\cdots+x_n^2\) の場合、
　　期待値 \(E(\chi^2(n))=n\) 、分散 \(V(\chi^2(n))=2n\) と予想されます。

一般の正規分布 \(N(\mu,\sigma^2)\)の場合、標準化した平方和がカイ二乗分布に従うので、
　　　 \(\chi^2=\left(\dfrac{x_1-\mu}{\sigma}\right)^2+\left(\dfrac{x_2-\mu}{\sigma}\right)^2+\cdots+\left(\dfrac{x_n-\mu}{\sigma}\right)^2\)
　　　 \(\chi^2=\dfrac{1}{\sigma^2}\left[(x_1-\mu)^2+(x_2-\mu)^2+\cdots+(x_n-\mu)^2\right]\)
偏差平方和を母分散で割ったものがカイ二乗値となります。

1.1.1 練習

独立な正規分布n個の偏差平方和
\(x_1,x_2,\cdots,x_5\sim N(1,2^2)\) に対するカイ二乗値は、
　　 \(\chi^2=\dfrac{1}{4}\left[(x_1-1)^2+(x_2-1)^2+\cdots+(x_5-1)^2\right]\)
偏差平方和の期待値を乱数を用いて求めましょう。

X = c()
for (i in 1:10000) {
  x <- rnorm(5,1,2)
  s <- sum((x-1)^2)
  X <-c(X,s)
}
mean(X)

## [1] 20.0365

偏差平方の期待値は分散 \(2^2\) であり、独立な５個の偏差平方の和なので、
期待値は \(2^2\cdot 5=20\) 　となります。

hist(X/4,prob=T)
curve(dchisq(x,5),add=T)

独立な５個の偏差平方の和を、母分散 \(2^2\) で割ったデータのヒストグラムは自由度５のカイ二乗分布に一致しています。

一般の \(\ S^2(n)=x_1^2+x_2^2+\cdots+x_n^2\) の場合、
　　期待値 \(E(S^2(n))=n\sigma^2\) 、分散 \(V(S^2(n))=2n\sigma^2\) と予想されます。

一般に、
独立なn個の偏差平方の和を母分散で割ったものは、n個の標準正規分布の平方和になるので、
自由度nのカイ二乗分布に従い、期待値はn、分散は2nになると考えられます。

1.2 自由度

上の例では \(X\sim N(0,1)\) というモデルが与えられており、 X の分散が 1 と分かっています。
しかし、実際の分析の場面では母分散も母平均も分かっておらず、得られたデータから推測する必要があります。
　　　\(V=\dfrac{x_1^2+x_2^2+\cdots+x_n^2}{n}=\dfrac{\chi^2}{n}\)　　　という式は、
母分散が未知の場合でも、得られたデータから標本分散が計算可能であることを示唆しています。

上の例では、
10個の \(x_i\sim N(0,1)\) が独立に動くとき、
　10個の和のカイ二乗値　\(\chi^2=x_1^2+x_2^2+\cdots+x_{10}^2\) 　は期待値が10、

これを10で割れば分散の値 1 が得られます。

5個の \(x_1,x_2,\cdots,x_5\sim N(1,2^2)\) が独立に動くとき、
　5個の偏差平方和　\(S^2=(x_1-1)^2+(x_2-1)^2+\cdots+(x_{5}-1)^2\) 　は期待値が \(5\cdot 2^2\) 、

これを 5で割れば分散の値 4 が得られます。

1.3 母平均が未知のとき

得られたデータだけから、母平均と母分散を求める方法を考えます。
　　（母平均）＝（標本平均）であることはよく知られています。
　　　\(E(X)=E(\bar{X})\) であるので、 　　 　　\(\mu\fallingdotseq \bar{X}=\dfrac{x_1+x_2+\cdots+x_n}{n}\)

近似値として等しいのですが、

　　\(S^2=(x_1-\mu)^2+(x_2-\mu)^2+\cdots+(x_n-\mu)^2\) 　のところを

　　\(S^2=(x_1-\bar{X})^2+(x_2-\bar{X})^2+\cdots+(x_n-\bar{X})^2\)　で代用した場合、
どれくらいズレが生じるかシミュレーションで見てみましょう。

X = c()
for (i in 1:10000) {
  x <- rnorm(10,0,1)
  s <- sum((x-mean(x))^2)
  X <-c(X,s)
}
mean(X)

## [1] 8.990354

ほとんど同じ式でありながら、カイ二乗値の期待値は10ではなく 9 になっています。
理由は、\((x_1-\bar{X}),~(x_2-\bar{X}),\cdots,(x_n-\bar{X})\) は独立でないため、自由度は(n-1)になっているのです。
一般の説明は難しいので、２個の偏差平方和の自由度が１であることを説明します。
\(x_1,x_2\sim N(0,1)\) は独立とします。
　　\(S^2=(x_1-\bar{X})^2+(x_2-\bar{X})^2=\dfrac{(x_1-x_2)^2}{2}\)　　 \(\left(\dfrac{x_1-x_2}{\sqrt{2}}\sim N(0,1)\right)\)
カイ二乗値は１個の標準正規分布の平方で表されるので、自由度は１となります。 シミュレーションをしてみると、

X = c()
for (i in 1:10000) {
  x <- rnorm(2,0,1)
  s <- sum((x-mean(x))^2)
  X <-c(X,s)
}
mean(X)

## [1] 0.9983033

hist(X,prob=T)

参考 高校生_t検定のpageファイルを参照

　　　\((x_1-\bar{X})^2+(x_2-\bar{X})^2=\dfrac{(x_1-x_2)^2}{2}\)
　　　\(\begin{split} E\left(\dfrac{(x_1-x_2)^2}{2}\right)&=\dfrac{1}{2}E\left(\left[(x_1-\mu) -(x_2-\mu)\right]^2\right)\\ &=\dfrac{1}{2}\left(E\left[(x_1-\mu)^2\right]+E\left[(x_2-\mu)^2\right]+E\left[(x_1-\mu)(x_1-\mu)\right]\right)\\ &=\dfrac{1}{2}\left(\sigma^2+\sigma^2+0\right)=\sigma^2 \end{split}\)

まとめ
\(S^2=(x_1-\bar{X})^2+(x_2-\bar{X})^2+\cdots+(x_n-\bar{X})^2\) の期待値は \((n-1)\sigma^2\)

このことは、偏差平方和が (n-1)個の独立な確率変数の平方和で表されることを示唆しています。

当サイト 正規分布と関連分布 参照

1.3.1 練習

平均3、標準偏差2の正規乱数10個の偏差平方和を、母平均を使わず 標本平均を使って計算します。

X = c()
for (i in 1:10000) {
  x <- rnorm(10,3,2)
  s <- sum((x-mean(x))^2)
  X <-c(X,s)
}
mean(X)/9

## [1] 4.039517

自由度は9なので、偏差平方和の期待値を9で割ったものが母分散の予測値になります。

1.4 不偏分散

母集団が \(X\sim N(\mu,\sigma^2)\) に従い、母集団から得られた標本 \(x_1,x_2,\cdots,x_n\)が独立に \(N(\mu,\sigma^2)\) に従うとします。
母平均 \(\mu\) が未知であるという前提で、母分散 \(\sigma^2\) の推定値を求めるには、

　　　\(\displaystyle u^2=\dfrac{1}{n-1}\sum_{i=1}^{n}{(x_i-\overline{x})^2}\)　これを不偏分散と言います。

根拠は、\(E\left(\sum_{i=1}^{n}{(x_i-\overline{x})^2}\right)=(n-1)\sigma^2\) であるからです。

2 Z検定

ヒストグラムや確率密度関数（確率分布のグラフ）は起こりやすさを表しています。
山の高いところは起こりやすく、山の低いところは起こりにくいところです。

• 平均を中心に左右標準偏差 2つ分(正確には1.96)に95％が入ると考えるのが統計的推測、
  
• 標準偏差 2つ分(正確には1.96)より中心から離れることはあり得ないと考えるのが統計的検定です。

2.1 標準正規分布

平均０, 標準偏差１の正規確率変数の確率密度関数は、
　　　\(\phi \left( z\right) =\dfrac{1}{\sqrt{2\pi }}\exp \left( -\dfrac{z^{2}}{2}\right)\sim N(0,1)\)　　で表され標準正規分布と言います。

curve(dnorm(x, 0, 1),xlim=c(-4, 4),ylim=c(0, 0.4), col=1,xlab="z")
curve(dnorm(x, 0, 1), from = qnorm(0.95, 0, 1), to = 4, #棄却域
      type = "h", col = "red", add = TRUE)
curve(dnorm(x, 0, 1), from = -qnorm(0.95, 0, 1), to = -4, #棄却域
      type = "h", col = "red", add = TRUE)

\(P(-1.96<z<1.96)=0.95\) （95%）であることから、\(|z|>1.96\) は起こり得ないと考えます。
図の赤い範囲を、確率的に起こりにくい範囲と考え、統計量がこの範囲に入った場合、帰無仮説を棄却します。( 棄却域といいます)

2.2 帰無仮説

Z検定は、母分散 \(\sigma^2\) が既知のとき、母平均が\(\mu\)であるという帰無仮説を検定します。
データが1個のとき
　　　\(X\sim N(\mu,\sigma^2)\) のもとで、データが有意かどうか検定するには、得られたデータが平均からどのくらい離れているかが問題になります。
得られたデータよりも平均から遠くなる確率をP値といい、関数 pnorm を利用して 求めることができます。
つまり、帰無仮説によって仮定された分布のもとで、データの統計量が起こりやすいかどうかを判断するのが仮説検定です。

• 右側に離れるかどうかを判断するのが右側検定（右側端5%が棄却域）
  
• 左側に離れるかどうかを判断するのが左側検定（左側端5%が棄却域）
  
• 両側に離れるかどうかを判断するのが両側検定（両側端2.5%が棄却域）

右側検定と左側検定をあわせて片側検定といいます。
まず、右側検定から見ていきます。
\(X\sim N(0,1^2)\) における有意水準5% 右側検定の境界は、X=1.64
この値より右に離れれば、あり得ないことが起こったと考えます。

2.2.1 母分散が\(1^2\)のとき

帰無仮説　\(X\sim N(0,1^2)\) のとき、

• データ X=1 は有意ですか？

curve(dnorm(x, 0, 1),xlim=c(-4, 4),ylim=c(0, 0.4), col=1,xlab="z")
curve(dnorm(x, 0, 1), from = qnorm(0.95, 0, 1), to = 4, #棄却域
      type = "h", col = "pink", add = TRUE)
curve(dnorm(x, 0, 1), from = 1, to = 1, #棄却域
      type = "h", col = "blue", add = TRUE)

X=1 は十分起こりうる値であり、帰無仮説は棄却されません。有意ではありません。

• データ X=2 は有意ですか？

curve(dnorm(x, 0, 1),xlim=c(-4, 4),ylim=c(0, 0.4), col=1,xlab="z")
curve(dnorm(x, 0, 1), from = qnorm(0.95, 0, 1), to = 4, #棄却域
      type = "h", col = "pink", add = TRUE)
curve(dnorm(x, 0, 1), from = 2, to = 2, #棄却域
      type = "h", col = "red", add = TRUE)

X=2は、起こりにくい範囲に入っており、帰無仮説は棄却されます。
P値は次の関数式で与えられ、

pnorm(2, mean=0, sd=1, lower.tail = FALSE)

## [1] 0.02275013

P値は 0.02 < 0.5　有意となります。

2.2.2 母分散が\(\sigma^2\)のとき

標準化式 　\(z=\dfrac{x-\mu}{\sigma}\sim N(0,1)\)
z はマハラノビスの距離といい、 x が平均から標準偏差いくつ分離れているかを表します。

これにより、
データが１個のとき
　　\(z=\dfrac{x-\mu}{\sigma}\sim N(0,1)\) により平均値の検定や区間推定ができます。
　　両側検定では \(|z|>1.96\) のとき、あり得ないことが起こったと考え平均\(\mu\) を棄却します。
　　右側検定では \(z >1.65\) のとき、あり得ないことが起こったと考え平均\(\mu\) を棄却します。

データがn個のとき
　　\(\bar{X}=\dfrac{x_1+x_2+\cdots+x_n}{n}\sim N\left(\mu,\dfrac{\sigma^2}{n}\right)\) により\(\bar{X}\)を標準化します。

　　\(z=\dfrac{\bar{x}-\mu}{\sigma/\sqrt{n}}\sim N(0,1)\) により平均値の検定や区間推定ができます。

2.2.3 データを見抜く（１）

母集団：\(X\sim N(0,1^2)\)から抽出された4個のサンプルが提出されました。
この標本が本当は\(X_i\sim N(1,1^2)~~~~(i=1:4)\) 平均１の母集団からのサンプルであるとして、そのことを見抜くことができるでしょうか。
検定を行ってみます。

統計量は　\(\bar{X}=\dfrac{x_1+x_2+\cdots+x_4}{4}\)

帰無仮説による分布は、\(\bar{X}\sim N\left(0,\dfrac{1}{4}\right)\)
サンプルの標本平均が棄却域に入るかどうかシミュレーションします。

x <- rnorm(4,1,1)
mean(x)

## [1] 1.306431

curve(dnorm(x, 0, 1/2),xlim=c(-4, 4),ylim=c(0, 0.8), col=1,xlab="z")#帰無仮説の分布
curve(dnorm(x, 0, 1/2), from = qnorm(0.95, 0, 1/2), to = 3, #棄却域
      type = "h", col = "pink", add = TRUE)
curve(dnorm(x, 0, 1/2), from = mean(x), to = mean(x), #棄却域
      type = "h", col = "blue", add = TRUE)

P値は

pnorm(mean(x), mean=0, sd=1/sqrt(4), lower.tail = FALSE)

## [1] 0.004489368

提出された偽データは棄却域に入り有意と判定されたでしょうか？
この判定を10000回繰り返し、判断の確度を見てみましょう。

2.2.4 検出力

X = c()
for (i in 1:10000) {
  x <- rnorm(4,1,1)
  p <- pnorm(mean(x), mean=0, sd=1/sqrt(4), lower.tail = FALSE)
  X <-c(X,p)
}
sum(X < 0.05)/10000

## [1] 0.6373

検出力は、対立仮説のもとで帰無仮説が棄却され有意となる確率です。
上ではシミュレーションによって検出力を求めましたが、数理的に計算することもできます。
対立仮説　\(X\sim N(1,0.5^2)\) 　のもとで

帰無仮説　\(X\sim N(0,0.5^2)\) 　が棄却される確率を求めます。

pnorm(1.645*0.5, mean=1, sd=1/2, lower.tail = FALSE)

## [1] 0.6387052

シミュレーションの結果と理論値は、ほぼ一致しています。

サンプルの数を100個に増やしてみます。

2.2.5 データを見抜く（２）

x <- rnorm(100,0.3,1)
mean(x)

## [1] 0.3373808

curve(dnorm(x, 0, 0.1),xlim=c(-2, 2),ylim=c(0, 4), col=1,xlab="z")#帰無仮説の分布
curve(dnorm(x, 0, 0.1), from = qnorm(0.95, 0, 0.1), to = 3, #棄却域
      type = "h", col = "pink", add = TRUE)
curve(dnorm(x, 0, 0.1), from = mean(x), to = mean(x), #棄却域
      type = "h", col = "blue", add = TRUE)

P値は

pnorm(mean(x), mean=0, sd=0.1, lower.tail = FALSE)

## [1] 0.0003706807

検出力は

pnorm(1.645*0.1, mean=0.3, sd=0.1, lower.tail = FALSE)

## [1] 0.9122912

このように標本数が増えると、わずかの差でも検出できるようになります。

3 t検定

t検定とは標準偏差に推定値（標準誤差）を用いる方法です。
t検定にはt分布が用いられます。
t分布はstudent分布と呼ばれ、ギネス醸造所の研究員であったゴセットがstudentというペンネームで研究を発表したことに由来します。 （ギネスビールと統計家参照）

数理的に次のことが成立します。（証明には立ち入りません）

\(W\sim \chi^2(m)\) 自由度mのカイ二乗分布に従い、標準正規分布に従う確率変数 Zが Wと独立ならば、
　　　\(t=\dfrac{Z}{\sqrt{\left. W \right/ m}} \sim t(m)\) 　　自由度mのt分布に従います。

t分布はどのような分布なのでしょうか。乱数によるシミュレーションによって調べてみましょう。
このテーマが高度であるのは、Zと Wが 互いに独立な確率変数、
従って、 Zと \(\sqrt{\left. W \right/ m}\) も 互いに独立な確率変数であり、t値は独立な確率変数の商の確率変数であるからです。

3.1 t分布

自由度10のカイ二乗乱数Wはrchisq(1, 10) によって 1個出力されます。
正規乱数をzとし、\(t=\dfrac{Z}{\sqrt{\left. W \right/ 10}}\)
によって t値を1000回出力させ、ヒストグラムを描きます。

X = c()
for (i in 1:10000) {
  t <- rnorm(1,0,1)/sqrt(rchisq(1, 10)/10)
  #t
  X <-c(X,t)
}
hist(X)

ZとWが互いに独立であることを確かめるため正規乱数とカイ二乗乱数の相関係数を計算してみます。

cor(rnorm(1000,0,1),rchisq(1000, 10))

## [1] -0.02786921

驚いたことに正規分布とよく似た形状になり、正規乱数とカイ二乗乱数は独立なので相関係数はほぼ０です。

標準正規分布の10個の平方和は自由度10のカイ二乗分布に従うので、この平方和を用いてシミュレーションしてみます。
縦軸を確率密度とし、標準正規分布の曲線(青) 、理論上のt分布曲線(赤) を比べてみます。

X = c()
for (i in 1:10000) {
  w <- rnorm(10,0,1)
  w <- sum(w^2)
  #w/10 #分散の近似値
  t <- rnorm(1,0,1)/sqrt(w/10)
  #t
  X <-c(X,t)
}
hist(X,prob=T)
curve(dnorm(x, 0, 1),xlim=c(-4, 4),ylim=c(0, 0.4), col="blue",add=T)
curve(dt(x, 10),xlim=c(-4, 4),ylim=c(0, 0.4), col="red",add=T)

t分布の確率密度関数は自由度を\(\nu\) として
　　　\(f(x) = \cfrac{\Gamma(\frac{\nu+1}{2})}{\sqrt{\nu\pi}\Gamma(\frac{\nu}{2})}\left(1+\cfrac{x^2}{\nu}\right)^{-\frac{\nu+1}{2}}\) 　　（\(\Gamma\) はガンマ関数）

で与えられますが、標準正規分布よりもガンマ関数のほうが当てはまりが良いようです。

また、自由度が大きくなると標準正規分布に近づいていきます。
下に、自由度100の場合のt分布を色を付けて示します。

X = c()
for (i in 1:10000) {
  t <- rnorm(1,0,1)/sqrt(rchisq(1, 10, ncp=0)/10)
  #t
  X <-c(X,t)
}

Y = c()
for (i in 1:10000) {
  t <- rnorm(1,0,1)/sqrt(rchisq(1, 100, ncp=0)/100)
  #t
  Y <-c(Y,t)
}

hist(X,prob=T,ylim=c(0,0.4))
hist(Y,prob=T,add=T,col="lightblue")
curve(dnorm(x,0,1),add=TRUE,col="red")

3.2 ゴセットの奇跡

ゴセッとの時代、コンピューターで乱数を発生させることはできませんでした。
佐藤達哉 様 によると、ゴセットは自宅で、特定の分布から小標本を繰り返し抽出し数値を計算することを手作業で繰り返したそうです。

標本平均を用いて計算したn個の偏差平方和 \(S^2\) を母分散\(\sigma^2\) で割ったUは、自由度(n-1)のカイ二乗分布に従うので、
　　　\(t=\dfrac{\overline{X}-\mu}{\left. s \right/ \sqrt{n}} =\dfrac{\left. (\overline{X}-\mu) \right/ \dfrac{\sigma}{\sqrt{n}}}{\left. s \right/ \sigma} =\dfrac{Z}{\sqrt{\left. s^2 \right/ \sigma^2}} =\dfrac{Z}{\sqrt{\left. S^2 \right/ (\sigma^2(n-1))}} =\dfrac{Z}{\sqrt{\left. U \right/ (n-1)}} \sim t(n-1)\)
（s は標本標準偏差）

標本平均 \(\overline{X}\) の標準化を標本分散で代用したものは、t分布に従うことが分かります。
（奇跡と呼ぶ理由）
母分散が未知の状況でデータのみが与えられているとき、母平均は \(\bar{X}\) によって 推定できますが、母分散が分からなければZ検定は使えず、母平均を検定・推定することができません。
標本偏差平方和は母平均が未知でも計算可能であり、帰無仮説に基づいた t値を算出すれば、そこから母平均についての検定をすることが可能になるということです。
偏差平方和 SとZが独立になるというのも驚くべきことです。
以下、n=4 の標本抽出を10000回行うシミュレーションを実施してみます。

Z = c()
W = c()
for (i in 1:1000) {
  x <- rnorm(4,1,1)
  s <- var(x)*3
  z <- (mean(x)-1)/sqrt(var(x))*2
  Z <-c(Z,z)
  W <-c(W,s)
}
hist(W)

hist(Z,prob=T,xlim=c(-4, 4),ylim=c(0, 0.4),breaks = "Scott")
curve(dt(x, 3),xlim=c(-4, 4),ylim=c(0, 0.4), col="red",add=T)
curve(dnorm(x, 0,1),xlim=c(-4, 4),ylim=c(0, 0.4), col="blue",add=T)

赤はt分布、青は正規分布ですが、あまり当てはまりはよくないようです。
これはデータ数が少ないためと考えられます。

データ数を100に増やしてみます。

  Z = c()
  W = c()
  for (i in 1:1000) {
    x <- rnorm(100,1,1)
    s <- var(x)*(99)
    z <- (mean(x)-1)/sqrt(var(x))*10
    Z <-c(Z,z)
    W <-c(W,s)
  }
  hist(Z,prob=T,xlim=c(-3, 3),ylim=c(0, 0.4),breaks = "Scott")
  curve(dt(x, 99),xlim=c(-4, 4),ylim=c(0, 0.4), col="red",add=T)
  curve(dnorm(x,0,1),xlim=c(-4, 4),ylim=c(0, 0.4), col="blue",add=T)

  mean(Z)

## [1] 0.003821635

  mean(W)

## [1] 99.3036

  cor(Z,W)

## [1] -0.006947886

かなり当てはまりは良くなりました。
検定が有効であるために、どのくらいのサンプル数が必要かは重大な問題です。
参考サイトを示しますが、今のところ、シミュレーションで当てはまりを確かめると良いでしょう。
t検定のサンプル数

3.3 サンプルサイズ

標本数を25にしてみます。

  Z = c()
  W = c()
  for (i in 1:1000) {
    x <- rnorm(25,1,1)
    s <- var(x)*(24)
    z <- (mean(x)-1)/sqrt(var(x))*5
    Z <-c(Z,z)
    W <-c(W,s)
  }
  hist(Z,prob=T,xlim=c(-3, 3),ylim=c(0, 0.4))
  curve(dt(x, 24),xlim=c(-4, 4),ylim=c(0, 0.4), col="red",add=T)
  curve(dnorm(x,0,1),xlim=c(-4, 4),ylim=c(0, 0.4), col="blue",add=T)

  mean(Z)

## [1] -0.01495919

  mean(W)

## [1] 24.0744

  cor(Z,W)

## [1] 0.01353488

ほぼ、当てはまっています。
（サンプル数を変えて確かめてみましょう）
サンプル数が10のとき

  Z = c()
  W = c()
  for (i in 1:1000) {
    x <- rnorm(10,1,1)
    s <- var(x)*(9)
    z <- (mean(x)-1)/sqrt(var(x))*sqrt(10)
    Z <-c(Z,z)
    W <-c(W,s)
  }
  hist(Z,prob=T,xlim=c(-3, 3),ylim=c(0, 0.4))
  curve(dt(x, 9),xlim=c(-4, 4),ylim=c(0, 0.4), col="red",add=T)
  curve(dnorm(x,0,1),xlim=c(-4, 4),ylim=c(0, 0.4), col="blue",add=T)

  mean(Z)

## [1] -0.05390614

  mean(W)

## [1] 8.897998

  cor(Z,W)

## [1] 0.02162459

3.4 1標本母平均検定（1群ｔ検定）

平均0.5のサンプル１０個を与えられたとします。
平均が０であるかどうかを、t検定によって検定します。

x <- rnorm(10,0.5,1)
t <- (mean(x)-0)/sqrt(var(x)/10)
pt(t, 9, lower.tail = FALSE)

## [1] 0.2071601

10000回繰り返して検出力を見ると、

X = c()
for (i in 1:10000) {
  x <- rnorm(10,0.5,1)
  t <- (mean(x)-0)/sqrt(var(x)/10)
  p <- pt(t, 9, lower.tail = FALSE)
  X <-c(X,p)
}
sum(X < 0.05)/10000

## [1] 0.4309

検出力は、あまり高くありません。
適切な判断のためには、サンプル数を多くする必要があります。

X = c()
for (i in 1:10000) {
  x <- rnorm(36,0.5,1)
  t <- (mean(x)-0)/sqrt(var(x)/36)
  p <- pt(t, 35, lower.tail = FALSE)
  X <-c(X,p)
}
sum(X < 0.05)/10000

## [1] 0.9064

サンプル数を36にすると検出力（検定が正しい確率）が90%を超えました。
一般に、１標本t検定ではサンプル数は30以上必要とされており、少ない場合はノンパラメトリックな検定法を用いるべきとされています。

２標本t検定、１標本t検定、ノンパラメトリックな検定法については、今後、応用ユニットで解説します。