1 カイ二乗値

確率変数 \(X\sim N(0,1)\) を自動販売機に喩えましょう。
ボタンを押すたびに値が出力され、その値はN(0,1)に従っています。

今あなたは小さな子供になって自販機の前に立っています。お母さんの言いつけで、自販機から値を持ってくるように言われたのです。あなたはその自販機から出る値の平均が０であることを、お母さんに聞いていました。でも、分散がどれだけかは知りません。
買い物かごを持ってきたあなたは、試しにボタンを押してみました。何回か押して、買い物カゴにたまった値を見てみると、０に近い値が多いのですが変動もあるようです。あなたはn回ボタンを押して、得られたn個の値から分散を求めることにしました。
母平均が０であることを予め知っているので、n個の出力値から計算した分散は、
　　　　　　\(V=\dfrac{x_1^2+x_2^2+\cdots+x_n^2}{n}=\dfrac{\chi^2}{n}\)

　　　（データ \((x_1,x_2,\cdots,x_n)\) は独立に標準正規分布に従う）

この分子の、 \(\chi^2=x_1^2+x_2^2+\cdots+x_n^2\)　　　をカイ二乗値といいます。

つまり、カイ二乗値は分散の計算式の分子ということになります。
（分散はカイ二乗値をnで割ったものと言うこともできます）
カイ二乗値は変動しうるもの、即ち確率変数であり、値を採取するたびにカイ二乗値は変動するということです。
読者は母分散が１であることを知っているので、カイ二乗値はnに近い値をとって変動するだろうと予測できます。
この変動がカイ二乗分布（自由度n）になります。

1.1 母平均が既知のとき

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.974443

var(X)

## [1] 19.97366

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

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

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

1.2 母平均が未知のとき

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

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

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

　　\(\chi^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] 9.040441

ほとんど同じ式でありながら、カイ二乗値の期待値は10ではなく 9 になっています。
理由は、\((x_1-\bar{X}),~(x_2-\bar{X}),\cdots,(x_n-\bar{X})\) は独立でないため、自由度は(n-1)になっているのです。

まとめ
\(\chi^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 不偏分散

母集団が \(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.313003

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.004319662

提出された偽データは棄却域に入り有意と判定されたでしょうか？
この判定を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.6314

検出力は、対立仮説のもとで帰無仮説が棄却され有意となる確率です。
上ではシミュレーションによって検出力を求めましたが、数理的に計算することもできます。
対立仮説　\(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.2775779

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.002753482

検出力は

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.03571753

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

標準正規分布の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))
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))
  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.004388903

  mean(W)

## [1] 98.71756

  cor(Z,W)

## [1] -0.008121231

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

  mean(W)

## [1] 23.89782

  cor(Z,W)

## [1] -0.0453193

ほぼ、当てはまっています。
（サンプル数を変えて確かめてみましょう）
サンプル数が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.01516656

  mean(W)

## [1] 8.953842

  cor(Z,W)

## [1] 0.03088135

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.07965904

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.429

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

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.9043

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

カイ二乗値とZ検定・t検定