1 不偏性と一致性.

分析では、データ \(X_1,\cdots,X_n\) から 母平均や母分散を推定することが必要になります。
また、その推定量が妥当であるかの検証も必要になります。
推定量の満たすべき妥当性は、

  • 不偏性：推定量の期待値が真値に一致する.
    （無限回シミュレーションした値の平均が真値と等しい）

  • 一致性：データ数が大きくなると、推定量が真値に近づく.
    （真値を中心としたヒストグラムの幅が 0 に近くなる）

「近づく」ことを数学では「収束する」といいますが、収束には様々な種類があり、ここでは収束の概念には立ち入らず直感的な理解にとどめます。

2 母平均の推定.

n個のデータが独立に得られているが、母平均と母分散は未知であるとしましょう。

まず、データ \(X_1,\cdots,X_n\) から 母平均を推定します。

2.1 相加平均.

標本平均 \(\bar{X}=\dfrac{X_1+X_2+\cdots +X_n}{n}\) は母平均の妥当な推定値と考えられます。

母集団は平均 0, 標準偏差 1 の正規分布に従うとし、データを rnorm(n,0,1) で出力します。

  • ４個のデータ rnorm(4,0,1) の標本平均を10000回出力したデータを M4.
  • 25個のデータ rnorm(25,0,1) の標本平均を10000回出力したデータを M25.

とします。

  • ヒストグラムの中心が、母平均 0 に近ければ不偏性が満たされ.
  • ヒストグラムの幅が、M4よりM25が小さくなっていれば一致性が満たされる.

と判断できます。
さらに、
標本平均の期待値、標本平均の標準偏差を mean関数 と sd関数 で計算します。

M4 = c()
for (i in 1:10000) {
  x <- rnorm(4,0,1)
  M4 <- c(M4,mean(x))
}

M25 = c()
for (i in 1:10000) {
  x <- rnorm(25,0,1)
  M25 <- c(M25,mean(x))
}
par(mfrow=c(1,2))
hist(M4,prob=T,xlim=c(-2,2))
hist(M25,prob=T,xlim=c(-2,2))

mean <- c(mean(M4),mean(M25))
sd <- c(sd(M4),sd(M25))
mean
## [1] -0.002684976  0.002199731
sd
## [1] 0.4968764 0.1983018

不偏性と一致性が満たされています.

2.2 相乗平均.

データが正である場合は、相乗平均を用いることも考えられます。

　　　\(\dfrac{x_1+\cdots+x_n}{n}\geq \sqrt[n]{x_1\cdots x_n}\).

であるので、相乗平均が小さめになることも予想されます。
確かめてみましょう。

母集団は、平均 5, 標準偏差 1 の正規分布に従うとし、データを rnorm(n,5,1) で出力します。

M4 = c()
for (i in 1:10000) {
  x <- rnorm(4,5,1)
  M4 <- c(M4,prod(x)^(1/length(x)))
}

M25 = c()
for (i in 1:10000) {
  x <- rnorm(25,5,1)
  M25 <- c(M25,prod(x)^(1/length(x)))
}
par(mfrow=c(1,2))
hist(M4,prob=T,xlim=c(3,7))
hist(M25,prob=T,xlim=c(3,7))

mean <- c(mean(M4),mean(M25))
sd <- c(sd(M4),sd(M25))
mean
## [1] 4.916631 4.899508
sd
## [1] 0.5136385 0.2072170

ヒストグラムは、相加平均に比べて少し左に寄っていないでしょうか？

「相乗平均＜相加平均」であるため、相乗平均の値は少し小さめに出力されているようです。

母平均は相加平均で計算されているので当然とは言えますが、
相乗平均は不偏性を満たさないので、推定量として適切ではありません。

3 母分散の推定.

データ \(X_1,\cdots,X_n\) から 母分散を推定します。

\(V(X)=E((X-\mu)^2)\) なので、分散は偏差平方の平均で推定できます。
しかし、その計算には母平均 μ が必要です。

  • 母平均 μ が既知の場合は、μ との偏差平方の平均で母分散を推定.

  • 母平均 μ が未知の場合は、\(\bar{X}\)との偏差平方の平均で母分散を推定.

このようになります。
上記の各方法で母分散を推定した場合、不偏性と一致性が満たされているかシミュレーションで確かめましょう。

3.1 母平均が既知のとき.

母集団は、平均 5, 標準偏差 1 の正規分布に従うとし、データを rnorm(n,5,1) で出力します。

  • ４個のデータ rnorm(4,5,1) の偏差平方の平均を1000回出力したデータを M4.
  • 25個のデータ rnorm(25,5,1) の偏差平方の平均を1000回出力したデータを M25.

とします。

偏差平方平均の期待値、偏差平方平均の標準偏差を mean関数 と sd関数 で計算します。

M4 = c()
for (i in 1:1000) {
  x <- rnorm(4,5,1)
  M4 <- c(M4,mean((x-5)^2))
}

M25 = c()
for (i in 1:1000) {
  x <- rnorm(25,5,1)
  M25 <- c(M25,mean((x-5)^2))
}
par(mfrow=c(1,2))
hist(M4,prob=T,xlim=c(0,4))
hist(M25,prob=T,xlim=c(0,4))

mean <- c(mean(M4),mean(M25))
sd <- c(sd(M4),sd(M25))
mean
## [1] 0.9626943 0.9744057
sd
## [1] 0.6829408 0.2733068

ヒストグラムは左右対称とは言えませんが、偏差平方平均の期待値は母分散 1 に近く、不偏性と一致性は満たされているようです。

3.2 母平均が未知のとき.

母集団は、平均 5, 標準偏差 1 の正規分布に従うことが分かっていますが、母平均が未知という前提なので、標本平均との偏差平方平均を計算します。

M4 = c()
for (i in 1:1000) {
  x <- rnorm(4,5,1)
  M4 <- c(M4,mean((x-mean(x))^2))
}

M25 = c()
for (i in 1:1000) {
  x <- rnorm(25,5,1)
  M25 <- c(M25,mean((x-mean(x))^2))
}
par(mfrow=c(1,2))
hist(M4,prob=T,xlim=c(0,4))
hist(M25,prob=T,xlim=c(0,4))

mean <- c(mean(M4),mean(M25))
sd <- c(sd(M4),sd(M25))
mean
## [1] 0.7504023 0.9531762
sd
## [1] 0.6023305 0.2664562

偏差平方平均の期待値は、母分散 1 より小さめに出ており、特に、データ数が４の M4 ではズレが大きいようです。
つまり、
母分散の推定値としては、偏差平方平均は不偏性が満たされていないようです。

どのくらいのズレが生じているか詳しく調べます。
個数で割ると変動が小さく判断しにくいので、偏差平方平均ではなく偏差平方和を調べます。
プログラムコード は、mean ではなく、 sum を用いることになります。

3.2.1 偏差平方和.

M4 = c()
for (i in 1:1000) {
  x <- rnorm(4,5,1)
  M4 <- c(M4,sum((x-mean(x))^2))
}

M25 = c()
for (i in 1:1000) {
  x <- rnorm(25,5,1)
  M25 <- c(M25,sum((x-mean(x))^2))
}
par(mfrow=c(1,2))
hist(M4,prob=T)
hist(M25,prob=T)

mean <- c(mean(M4),mean(M25))
mean
## [1]  2.935337 23.980829

データ数との関連が見えてきました。
偏差平方和の期待値は（n−１）となっており、母分散が \(\sigma^2=1\) であることから.

　　　\((X_1-\bar{X}),(X_2-\bar{X}),\cdots ,(X_n-\bar{X})\) の平方和の期待値は \((n-1)\sigma^2\).

と考えられます。.

3.3 課題.

分散を変えて、このことを確かめましょう。

3.3.1 平均５, 分散４.

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

M25 = c()
for (i in 1:10000) {
  x <- rnorm(25,5,2)
  M25 <- c(M25,sum((x-mean(x))^2))
}
par(mfrow=c(1,2))
hist(M4,prob=T)
hist(M25,prob=T)

mean <- c(mean(M4),mean(M25))
mean
## [1] 12.05492 95.81887

　　　\(12=(4-1)\times 2^2\)　　　　\(96=(25-1)\times 2^2\).

偏差平方和の期待値は 　\((n-1)\sigma^2\)　と考えられます。.
即ち、
　　　\((X_1-\bar{X})^2+(X_2-\bar{X})^2+\cdots +(X_n-\bar{X})^2\) 　の期待値は \((n-1)\sigma^2\).

さらに、
　　　\(U=\dfrac{(X_1-\bar{X})^2+(X_2-\bar{X})^2+\cdots +(X_n-\bar{X})^2}{\sigma^2}\) 　　とおくと、

その期待値は \((n-1)\)　です。

3.3.2 研究課題.

数学的証明に関心のある人は、２個の偏差平方和について、数理的導出を考えてみましょう。

\(X_1,X_2\sim N(0,1)\) は独立とします。
\((X_1-\bar{X}),(X_2-\bar{X})\) の平方和の期待値は \((2-1)\cdot 1^2=1\)と予想されます。.

まず、シミュレーションで確かめます。

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

hist(M2,prob=T)

mean <- c(mean(M2))
mean
## [1] 0.9873602

次に、式の変形で説明します。

\(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)\)
2 個の偏差平方和は１個の標準正規分布の平方で表され、その期待値は１となります。

個数が n 個になった場合も、n 個の偏差平方和は(n-1)個の標準正規分布の平方で表され、その期待値は(n-1)となります。
この導出にはヘルマート行列による一次変換が必要です。
線形代数で「行列と一次変換」を学んでから挑戦してください。
（ひょっこり統計参照）.

証明には立ち入りませんが、

　　　\(U=\dfrac{(X_1-\bar{X})^2+(X_2-\bar{X})^2+\cdots +(X_n-\bar{X})^2}{\sigma^2}\)

は、(n-1)個の独立な標準正規分布の平方和に変形可能であり、自由度 (n-1) のカイ二乗分布に従うことになります。

3.3.3 まとめ.

偏差平方和の期待値は、母分散の (n-1)倍であることが分かりました。
(n-1) を偏差平方和の自由度といいます。

偏差平方和の自由度は、

  • 数理的には、何個の独立な確率変数の平方和で表されるか、その個数ということができます。
  • シミュレーションでは、偏差平方和の期待値が母分散の何倍であるか、その倍数ということができます。

偏差平方和を (n-1) で割れば、その期待値は母分散になることを確かめておきます。

M4 = c()
for (i in 1:10000) {
  x <- rnorm(4,5,2)
  M4 <- c(M4,sum((x-mean(x))^2)/3)
}

M25 = c()
for (i in 1:10000) {
  x <- rnorm(25,5,2)
  M25 <- c(M25,sum((x-mean(x))^2)/24)
}
par(mfrow=c(1,2))
hist(M4,prob=T)
hist(M25,prob=T)

mean <- c(mean(M4),mean(M25))
sd <- c(sd(M4),sd(M25))
mean
## [1] 3.929045 4.020937
sd
## [1] 3.197256 1.174435

確かに、(n-1)で割れば偏差平方和の期待値は母分散になりますが、データ数が小さいときは sd 、即ち、誤差が大きすぎるようです。

この推定値を用いて標準化を行なっても標準正規分布にはなりません。
次の章で扱う t分布 になります。

4 t分布.

t分布はstudent分布と呼ばれ、ギネス醸造所の研究員であったゴセットがstudentというペンネームで研究を発表したことに由来します。 （ギネスビールと統計家参照）

製造業では製品の検査をして出荷しますが、検査による製品ロスを少なくするために、検査する標本数は小さくする必要があります。
しかし、標本数が小さければ、標準化したデータの標準正規分布とのズレが大きくなります。小標本の場合に、どのようなズレが生じるのか調べておく必要があります。
ゴセットの時代、コンピューターで乱数を発生させることはできませんでした。
佐藤達哉 様 によると、ゴセットは自宅で、特定の分布から小標本を繰り返し抽出し数値を計算することを手作業で繰り返したそうです。
ゴセットは、
　　　「標本平均 \(\bar{X}\) の標準化を標本分散で代用したものは、t分布に従う」

ということを発見します。
詳しくは次の章で見ていきましょう。