3章 母平均の推定

独立な観測データ \(X_1,\cdots,X_n\) が得られたとき、母平均 μ を区間推定することを目的とします。その際、

• 母分散が既知のとき.
  
• 母分散が未知のとき.

を別々に考える必要があります。

まず、母分散が既知のときから考えます。

1 母分散が既知.

1.1 標準誤差.

独立な観測データ \(X_1,\cdots,X_n\) が得られたとすると、
　　　その標本平均 \(\bar{X}\) を推定値と考えるのは自然なことです。
同時に、推定精度がどのくらいかも把握しておかなければなりません。

確率モデルとして \(X_1,\cdots,X_n \sim N(\mu,\sigma^2)\) は独立に得られたデータとすると、
　　　 \(\bar{X}=\dfrac{X_1+X_2+\cdots +X_n}{n}\sim N\left(\mu,\dfrac{\sigma^2}{n}\right)\).
２章 「2.1 相加平均」は、この公式が成立する具体例となっています。

• 各観測データ \(X_i\) のばらつきを、標準偏差 \(\sigma\)
  
• 標本平均 \(\bar{X}\) のばらつきを、標準誤差 \(\dfrac{\sigma}{\sqrt{n}}\)

と区別して呼びます。母平均の推定値の標準偏差が標準誤差ということになります。
この標準誤差が推定精度です。

母分散が既知であれば、

• \(\dfrac{\bar{X}-\mu}{\left. \sigma \right/ \sqrt{n}}\sim N(0,1)\)　　が成立します。
  標本平均を標準誤差で標準化したものです。

この式から、観測データ \(X_1,\cdots,X_n\)から得られる \(\bar{X}\) を用いて、母平均 μ を区間推定することが可能になります。即ち、
95％の確率で 　\(-1.96<\dfrac{\bar{X}-\mu}{\left. \sigma \right/ \sqrt{n}}<1.96\) 　であることから、

• 95％の確率で　\(\bar{X}-1.96\cdot\dfrac{\sigma}{\sqrt{n}}<\mu<\bar{X}+1.96\cdot \dfrac{\sigma}{\sqrt{n}}\).

これは、標本平均 \(\bar{X}\) の実現値を中心として左右に標準誤差の1.96倍幅の区間に母平均 μ が存在することを意味し、母平均 μ の95％信頼区間といいます。

1.2 信頼区間.

信頼区間は、標本平均の実現値と標準誤差によって表示できます。

　　　dnorm(x, 標本平均の実現値, 標準誤差).

をcurve( ) で描画したとき、実現値を中心として左右に標準誤差の1.96倍幅の区間が信頼区間です。ざっくり言えば、山になっている部分が信頼区間といってもよいでしょう。

「25個のデータの標本平均 mean(rnorm(25,0,1)) をもとに μ の信頼区間を表示する」

curve( )関数を用いて、信頼区間を3回表示します。

par(mfrow=c(3,1))
curve(dnorm(x, mean(rnorm(25,0,1)), 1/5),xlim=c(-2,2),type = "h", col = "red",ylab="")
curve(dnorm(x, mean(rnorm(25,0,1)), 1/5),xlim=c(-2,2),      type = "h", col = "red",ylab="")
curve(dnorm(x, mean(rnorm(25,0,1)), 1/5),xlim=c(-2,2),      type = "h", col = "red",ylab="")

［参考］

正確に95%信頼区間を表示するには、コードに「from to」を加えますが、煩雑になります。上記のコードで十分イメージできると思います。

par(mfrow=c(3,1))
m <- mean(rnorm(25,0,1))
curve(dnorm(x, m, 1/5),xlim=c(-2,2),type = "h", col = "red",ylab="",from = qnorm(0.025,m, 1/5), to = qnorm(0.975, m, 1/5))
m <- mean(rnorm(25,0,1))
curve(dnorm(x, m, 1/5),xlim=c(-2,2),type = "h", col = "red",ylab="",from = qnorm(0.025,m, 1/5), to = qnorm(0.975, m, 1/5))
m <- mean(rnorm(25,0,1))
curve(dnorm(x, m, 1/5),xlim=c(-2,2),type = "h", col = "red",ylab="",from = qnorm(0.025,m, 1/5), to = qnorm(0.975, m, 1/5))

赤の領域に対応する x の範囲が推定区間と考えることができます。

• 得られた標本平均の実現値に応じ推定区間が移動する.
  
• 推定区間は母平均 0 を含んでいる（正確には、95%の確率で）.

ということが読み取れます。

「この操作を100回実施すると、95回は母平均 0 を含む」.
これが95%信頼区間の意味です。

95%信頼区間には母平均を含まないリスクが 5% あることになり、リスクを意識した表現といえます。
また、推定区間内の確率密度は均一ではないので、このように確率密度を意識した表現もあっていいでしょう。

2 母分散が未知.

母分散が未知のときは、どうしたらよいのでしょうか？
２章より、母分散の推定値として

　　　\(s^2=\dfrac{(X_1-\bar{X})^2+(X_2-\bar{X})^2+\cdots +(X_n-\bar{X})^2}{n-1}\).
は不偏性と一致性を満たしており、不偏分散と呼ばれます。
しかし、\(s^2\)は、データ\(X_1,\cdots,X_n\)を得るごとに変動する確率変数です。
つまり、

• \(\dfrac{\bar{X}-\mu}{\left. s \right/ \sqrt{n}}\) は確率変数の商.

であるので、どのような分布になるかは分かりません。
調べて見ましょう。

この実験は、観測データ \(X_1,\cdots,X_n\) の実現値から母平均が推定可能かどうかを調べるものです。
データには、母数（母平均、母分散）が既知のデータを生成して用いますが、
推定には n 個の実現値のみを用い、母分散は使いません。

• 母分散が既知のときは.
  \(\dfrac{\bar{X}-\mu}{\left. \sigma \right/ \sqrt{n}}\) が標準正規分布に従うことから区間推定が可能でした。

• 母分散が未知のとき.
  \(\dfrac{\bar{X}-\mu}{\left. s \right/ \sqrt{n}}\) が従う分布を調べ、区間推定に利用します。
  （これはゴセットが発見した分布で、t分布といいます）

　\(\dfrac{\bar{X}-\mu}{\left. s \right/ \sqrt{n}}\) を本書では ゴセットの標準化、分布をゴセットのt分布と呼ぶことにします。

2.1 ゴセットの t分布.

母集団は、平均 5, 標準偏差 2 の正規分布に従うとします。
データ rnorm(n,5,2) に対するゴセットの標準化は、

• \(\dfrac{\bar{X}-\mu}{\left. s \right/ \sqrt{n}}=\left(\dfrac{mean(x)-5}{sd(x)}\right)\cdot \sqrt{n}\) で計算されます。

分母が sd(x) という関数のままであることに注意してください。

• ４個のデータ x=rnorm(4,5,2) に対し、ゴセットの標準化を1000回出力したデータを S4.
  
• 25個のデータ rnorm(25,5,2) に対し、ゴセットの標準化をを1000回出力したデータを S25.

とします。

また、

• 標準正規分布 dnorm(0,1) を赤の曲線で、
  
• ゴセットが発見した t分布の曲線 dt(x, 標本数−１) を青で、

描画します。

S4 = c()
for (i in 1:1000) {
  x <- rnorm(4,5,2)
  S4 <- c(S4,(mean(x)-5)/sd(x)*2)
}

S25 = c()
for (i in 1:1000) {
  x <- rnorm(25,5,2)
  S25 <- c(S25,(mean(x)-5)/sd(x)*5)
}
par(mfrow=c(1,2))
hist(S4,prob=T,xlim=c(-5,5),ylim=c(0,0.4),breaks="Scott")
curve(dnorm(x, 0, 1),add=T,col="red")
curve(dt(x, 3),add=T,col="blue")

hist(S25,prob=T,xlim=c(-5,5))
curve(dnorm(x, 0, 1),add=T,col="red")
curve(dt(x, 24),add=T,col="blue")

標本数が４のときは、t分布のほうが当てはまりが良いことが分かります。
t分布の裾が高いということは、推定区間の幅を広げる必要があることを示唆します。
同じ推定を使うと、確率が95%にならないということです。
少数標本で標準正規分布を用いると、推定のリスクが変わるので適切ではありません。
この理由で、少数標本ではゴセットのt分布を用いた推定が用いられています。

2.2 信頼性.

t分布を使わず標準正規分布で推定した場合、信頼性はどの程度下がるでしょうか？

標準正規分布における95%信頼区間の境界値は 1.96 です。
これを t分布に当てはめ、信頼性がどの程度下がるか見てみましょう。

標本数が４のときの信頼性 P4、標本数が 25 のときの信頼性 P25 を求めます。

P4 <- pt(1.96,3)-pt(-1.96,3)
P25 <- pt(1.96,24)-pt(-1.96,24)
c(P4,P25)

## [1] 0.8551478 0.9382940

標本数が４の場合、リスクが 15% あり、適当とはいえません。
標本数が25になると信頼性は95%に近く、リスクは小さいようです。

標本数を４から３０まで変化させたとき、信頼性がどう変わるかグラフにします。

P <- function(n) {
  p <- pt(1.96,n-1)-pt(-1.96,n-1) 
  return(p)
}
n <- 4:30
plot(n,P(n))

正規分布で近似するには、標本数は30以上必要とされています。
信頼度のグラフを見ると、標本数が30を越えれば標準正規分布で近似しても差し支えないようです。

これは、中心極限定理(Central Limit Theorem) から、 標本数が大きくなると標本平均の分布は正規分布に近づくことによります。

2.3 t分布の定義.

数学では次のような方法がよく用いられます。

• 一般化することで、物事を簡明化する.
  
• 簡明化されたことを、様々な分野に応用する.

ゴセットによって発見された t分布も、一般化され数理統計のテキストでは次のように定義されています。

［定義］

• 標準正規分布に従う \(Z\sim{}N(0,1)\).
• 自由度mのカイ二乗分布に従う \(U\sim{}\chi^2(m)\).

について、\(Z,~U\) が独立であるとき.
　　\(T=\dfrac{Z}{\sqrt{\left. U \right/ m}}\)　の従う分布を 自由度mのt分布という。

確率密度関数は、
　　\(f(x) = \cfrac{\Gamma(\frac{m+1}{2})}{\sqrt{m\pi}\Gamma(\frac{m}{2})}\left(1+\cfrac{x^2}{m}\right)^{-\frac{m+1}{2}}\) （\(\Gamma\) はガンマ関数）.

で与えられる。

数理的証明は難しいので、「t分布の定義」と「ゴセットの標準化」のつながりを確認するにとどめます。あまり深入りせず読んでください。

［カイ二乗値］.

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

即ち、
標本平均を用いて計算したn個の偏差平方和を母分散 \(\sigma^2\) で割った統計量は、

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

\(U\) は自由度 (n-1) のカイ二乗分布に従います。

母分散の推定量は、

　　　\(s^2=\dfrac{(X_1-\bar{X})^2+(X_2-\bar{X})^2+\cdots +(X_n-\bar{X})^2}{n-1}=\dfrac{\sigma^2 U}{n-1}\).
　　　\(t=\dfrac{\bar{X}-\mu}{\left. s \right/ \sqrt{n}} =\dfrac{\left. (\bar{X}-\mu) \right/ \dfrac{\sigma}{\sqrt{n}}}{\left. s \right/ \sigma} =\dfrac{Z}{\sqrt{\left. s^2 \right/ \sigma^2}} =\dfrac{Z}{\sqrt{\left. U \right/ (n-1)}} \sim t(n-1)\)
（s は標本標準偏差）

「標準正規分布 Z 」と「カイ二乗分布／自由度　の平方根」の商に変形され、ゴセットの標準化は 自由度(n-1)のt分布に従うことが分かります。

［まとめ］.

• ゴセットはゴセットの標準化を用いてゴセットのt分布を；発見しました。（自由度は n-1）.
• フィッシャーはゴセットの発見の重要性を見抜き、数学的に定義し証明しました。自由度nのカイ二乗分布を用いた定義ではt分布の自由度もnとなります）。
  
• ゴセットのt分布は自由度(n-1)のカイ二乗分布を用いるので、自由度は(n-1)になります。

t分布の数理は奥が深く難しい印象です。

ゴセットのt分布は自由度が (n-1) であることに注意し、一般のt分布と区別したほうが良いでしょう

（ゴセットとフィッシャー）
ゴセットは1908年、ビール酵母の計測誤差についての論文で、ゴセットの標準化によるt分布曲線を示しましたが、フィッシャーが提案した数理的証明には理解を示さなかったようです。数理的定義と現実との乖離をゴセット自身も感じていたようです。

2.4 信頼区間.

観測データ \(X_1,\cdots,X_n\) の実現値から求めたゴセットの標準化.
　　　\(t=\dfrac{\bar{X}-\mu}{\left. s \right/ \sqrt{n}}\) 　が t分布に従うとき.

95%信頼区間の境界値は、自由度(n-1)に従って変化します。
境界値を　a <- qt(0.975, n-1)　によって算出すると、

95%の確率で 　\(-a<\dfrac{\bar{X}-\mu}{\left. s \right/ \sqrt{n}}<a\) が成立するので.
　　　\(\bar{X}-a{\dfrac{s}{\sqrt{n}}}<\mu<\bar{X}＋a{\dfrac{s}{\sqrt{n}}}\).

これが、t分布による信頼区間です。

信頼区間の幅に影響する95%境界値について、
　「自由度 3 のときの境界値 a3 　自由度 24 のときの境界値 a24」
を求めてみると、

a3 <- qt(0.975, 3)
a24 <- qt(0.975, 24)
a <- c(a3,a24)
a

## [1] 3.182446 2.063899

自由度が 3 、即ち、標本数が 4 のときは、信頼区間の幅が広くなっています。

　　　\(t=\dfrac{\bar{X}-\mu}{\left. s \right/ \sqrt{n}}\) 　より　\(\bar{X}=\dfrac{s}{\sqrt{n}}\cdot t+\mu\)

\(\bar{X}\)の確率密度関数のグラフは t分布の関数 t(x,n-1) を用いて描画できます。
ここでは証明には立ち入らず、

• \(Z=\dfrac{X-\mu}{\sigma}\) の確率密度関数 　\(\phi \left( z\right) =\dfrac{1}{\sqrt{2\pi }}\exp \left( -\dfrac{z^{2}}{2}\right)\).
  
• \(X=\sigma Z+\mu\)の確率密度関数　\(f(x)=\dfrac{1}{\sqrt{2\pi\sigma}}\exp\left(-\dfrac{(x-\mu)^ 2}{2\sigma^ 2}\right)\)

から類推して用います。

t分布による信頼区間は、標本が得られるたびに標準誤差が変化するので、信頼区間の幅が変化します。

2.4.1 標本数４.

par(mfrow=c(3,1))
data1 <- rnorm(4,0,1)
curve((2/sd(data1))*dt((2/sd(data1))*(x-mean(data1)), 3),xlim=c(-2,2),
      ylim=c(0,1),type = "h", col = "blue",ylab="")
data2 <- rnorm(4,0,1)
curve((2/sd(data2))*dt((2/sd(data2))*(x-mean(data2)), 3),xlim=c(-2,2),
      ylim=c(0,1),type = "h", col = "blue",ylab="")
data3 <- rnorm(4,0,1)
curve((2/sd(data3))*dt((2/sd(data3))*(x-mean(data3)), 3),xlim=c(-2,2),
      ylim=c(0,1),type = "h", col = "blue",ylab="")

幅とその変化が大きいので、標本数４は小さすぎるという感じがします。

［偏差平方和 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)
cor(Z,W)

## [1] -0.01440749

分母と分子は無相関になっており、独立と考えられます。

2.4.2 標本数16.

par(mfrow=c(3,1))
curve(dt((4/sd(rnorm(16,0,1)))*(x-mean(rnorm(16,0,1))), 8),xlim=c(-2,2),
      ylim=c(0,0.5),type = "h", col = "blue",ylab="")
curve(dt((4/sd(rnorm(16,0,1)))*(x-mean(rnorm(16,0,1))), 8),xlim=c(-2,2),
      ylim=c(0,0.5),type = "h", col = "blue",ylab="")
curve(dt((4/sd(rnorm(16,0,1)))*(x-mean(rnorm(16,0,1))), 8),xlim=c(-2,2),
      ylim=c(0,0.5),type = "h", col = "blue",ylab="")

  Z = c()
  W = c()
  for (i in 1:1000) {
    x <- rnorm(16,1,1)
    s <- var(x)*(15)
    z <- (mean(x)-1)/sqrt(var(x))*sqrt(16)
    Z <-c(Z,z)
    W <-c(W,s)
  }
  cor(Z,W)

## [1] 0.009019971

標本数が増えると、区間幅の変動は小さくなり区間推定は妥当なものになってきます。

一般に、標本数は30以上と言われています。
30以上であれば、正規分布にも近くなり、精度のよいものと言えるでしょう。
少ない標本数しかない場合は、シミュレーションで変動をよく確かめる必要があるでしょう。