医学情報処理演習第４回「標本統計量と母数推定」 2006 年 10 月 30 日 中澤

港¶
前回の課題の回答例³sleep データの分布の位置とばらつきの情報を与える記述統計量を計算する前に，図示により分布の様子を確認する。図示の方法はいろいろあるが，とりあえずヒストグラムを作ってみる。

¶³attach(sleep)
layout(c(1,2))
hist(extra[group=="1"])hist(extra[group=="2"])detach(sleep)µ´

により下図ができ，正規分布とは到底言えなそうだとわかる。そこで位置とばらつきの情報としては，平均±不偏標準偏差よりもむしろ中央値±四分位偏差を示すべきと判断される。

¶³

SIQR <- function(x) { (fivenum(x)[4]-fivenum(x)[2])/2 }attach(sleep)
tapply(extra,group,median)tapply(extra,group,SIQR)
detach(sleep)µ´

により（tapply は，１番目の引数について，２番目の引数ごとに層別して，３番目の引数で与える関数を適用するという関数なので，group 別の統計量を計算する時に便利），催眠薬１を投与したときの睡眠時間変化の中央値±四分位偏差は 0.35 ± 1.1，催眠薬２投与時の睡眠時間変化の中央値±四分位偏差は1.75 ± 1.8 とわかる。
µ´

標本統計量と母数前回も軽く説明したように，我々が統計解析で相手にするのは標本データの場合が多く，得られる統計量も標本統計量であるが，本当に知りたいのは母集団の統計量（母数）であるため，標本統計量から母数を推定しなければならない。その仕組みを考えるため，まず標本抽出をシミュレートしてみる。

1

標本抽出
表 1. 成人男子の身長からなる母集団と母平均 µ，母分散 σ 2 の計算
身長 (cm) (x)
度数 相対度数（p(x))
x · p(x) (x − µ)2 · p(x)

144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
合計

2
3
6
17
19
56
125
219
463
915
1609
2649
4550
7214
11005
16081
22098
29903
39048
48312
57703
66639
73332
78051
79829
77866
73767
66321
57993
48410
39081
29967
22055
15810
10875
7309
4596
2726
1519
939
462
224
128
50
31
14
5
4
1000000

0.000002
0.000003
0.000006
0.000017
0.000019
0.000056
0.000125
0.000219
0.0004630.000915
0.001609
0.002649
0.004550
0.007214
0.011005
0.016081
0.022098
0.029903
0.039048
0.048312
0.057703
0.066639
0.073332
0.078051
0.079829
0.077866
0.073767
0.066321
0.057993
0.048410
0.039081
0.029967
0.022055
0.015810
0.010875
0.007309
0.004596
0.002726
0.001519
0.000939
0.000462
0.000224
0.0001280.000050
0.000031
0.000014
0.000005
0.000004
1.00

0.000288
0.000435
0.000876
0.002499
0.002812
0.008344
0.018750
0.033069
0.070376
0.139995
0.247786
0.410595
0.709800
1.132598
1.738790
2.556879
3.535680
4.814383
6.325776
7.874856
9.463292
10.995435
12.173112
13.034517
13.411272
13.159354
12.540390
11.340891
9.974796
8.374930
6.800094
5.244225
3.881680
2.798370
1.935750
1.308311
0.827280
0.493406
0.276458
0.171837
0.085008
0.041440
0.023808
0.009350
0.005828
0.002646
0.000950
0.000764µ = 167.9998

2

1.151979e−031.586970e−03
2.903942e−03
7.496844e−03
7.599834e−03
2.021553e−02
4.049901e−02
6.328937e−021.185248e−01
2.058690e−01
3.153541e−01
4.476659e−01
6.551761e−01
8.728592e−01
1.100452e+00
1.302498e+00
1.414195e+00
1.465155e+00
1.405625e+00
1.207694e+00
9.231469e−01
5.996634e−01
2.932638e−01
7.801682e−02
3.828679e−09
7.790011e−02
2.951326e−01
5.969761e−01
9.279896e−01
1.210356e+00
1.407019e+00
1.468475e+00
1.411597e+00
1.280672e+00
1.087548e+00
8.844242e−01
6.618482e−01
4.607095e−01
2.977333e−01
2.112812e−01
1.182752e−01
6.473767e−02
4.147301e−02
1.805042e−021.240027e−02
6.174129e−03
2.420048e−03
2.116040e−03
σ 2 = 25.09521
σ = 5.009512

100 万人の成人男性の身長からなる母集団が，表１のようになっていたとしよう∗1 。

表１のように計算すれば，X の平均 µ =Px · p(x) = 167.9998，X の分散 σ 2 = 25.09521 である。ここから 10 個の標本を抽出することを考える。なるべく母集団全体の情報を偏りなく代表させるためには，身長が書かれた 100 万枚のカードを袋に入れてよくかきまぜ，10 枚取り出すことが考えられる。無作為に抽出される最初のカードが x である確率は，p(x) である。この 1 回目の抽出を X1 であらわす。2 枚目のカードについても同じである（厳密に言えば，1 枚とっては元に戻して再びよくかきまぜて抽出する「復元抽出」でないと，1 枚目に何が出たかに依存して相対度数が変化するため，条件付き確率はまったく同じにはならないが，100 万くらいの大きな母集団からサイズ 10 の標本を取り出すのであれば，「非復元抽出」（1 度サンプルを取り出したら元に戻さない）であっても，p(x) はほとんど変化しないし，1 枚目に何がでるかという場合を総当りして合計すれば，2 枚目に何が出るかという確率も厳密に p(x) となる）。このように，無作為に選び出されるカードは，すべて母集団における相対度数によって与えられる確率で，母集団の値のどれかをとりうる，確率変数である。
1(X1 + X2 + ... + X10 ) なので，その期待値は，サイズ 10 の標本の標本平均 X̄ は，X̄ = 10E(X̄) =111(E(X1 ) + E(X2 ) + ... + E(X10 )) =(µ + µ + ...µ) =· 10 · µ = µ101010となって母平均に一致する。標本平均の分散はV (X̄) = V ({1(X1 + X2 + ...X10 )})101 2) {V (X1 ) + V (X2 ) + ... + V (X10 )}1011= ( )2 (σ 2 + σ 2 + ... + σ 2 ) = ( )2 · 10 · σ 21010=(=σ210したがって，標本平均の標準偏差は，σX̄ = √σ10 となる。これは，標本から推定される母集団の平均値の信頼性を示す値なので，平均の標準誤差と呼ばれる。標本の標準偏差が生データのばらつきを示す値であってサン∗1 こんなにきれいな正規分布に近い表が現実のデータとして存在するはずがなく，実は，これは R の正規乱数を使って得た人工データである。母平均 168，母標準偏差 5 の正規乱数を 100 万個発生させるためには，X <- rnorm(1000000,168,5) とすればよい。
実行するたびに異なるデータになってしまうのを避けるためには，このコマンドを実行する前に乱数発生アルゴリズムとその初期値を RNGkind("Mersenne-Twister",normal.kind="Inversion"); set.seed(1) などとして固定すればよい。

3

プルサイズに依存しないのに対し，標準誤差は平均値の信頼性を示す値なので，一般にサンプルサイズを大きくすれば小さくできる値であるため，標準誤差を示す時はサンプルサイズも明記しなくてはならない。
一方，母集団の分散の推定値としての標本の不偏分散は，
10V (X) =1X(Xi − X̄)29 i=11である。なぜ 10でなくて 19 かというと，X̄ は 10 個の Xi からの計算値なので，偏差の情報は見かけ上 10 個あるようにみえるが，9 個分しかないからである（10 個の値の平均と 9 個の値が決まれば 10 個目の値はその計算値として示せてしまう）。厳密な証明はやや高度なので省略するが（参考書としては，竹村彰通『現代数理統計学』創文社をお薦めする），以下，数値シミュレーションで試してみよう。

R で標本抽出をする関数は，sample() である。replace=オプションで復元，非復元を選択できる（デフォルトは非復元）。sample() 関数は乱数を使って標本抽出しているので，乱数を初期化して揃えない限り，実行するたびに結果は変わる。

¶it04-1-2006.R³

HT <- 144:191
NUM <- c(2,3,6,17,19,56,125,219,463,915,1609,2649,4550,7214,11005,16081,22098,29903,39048,48312,57703,66639,73332,78051,79829,77866,73767,66321,57993,48410,39081,29967,22055,15810,10875,7309,4596,2726,1519,939,462,
224,128,50,31,14,5,4)
X <- rep(HT,NUM)
RNGkind("Mersenne-Twister",normal.kind="Inversion")set.seed(1)layout(c(1,2,3,4))
hist(s10 <- sample(X,10))
hist(s100 <- sample(X,100))hist(s1000 <- sample(X,1000))
hist(s10000 <- sample(X,10000))
print(c(mean(s10),mean(s100),mean(s1000),mean(s10000)))
print(c(sd(s10),sd(s100),sd(s1000),sd(s10000)))µ´

サンプルサイズが大きくなるにつれて，標本分布は母集団の分布に近づく（平均値からみても，不偏標準偏差からみても）。もっとも，ある程度大きなサンプルになると，それほどの差は感じられない。ここで標本の4分散でなくて不偏分散の方が母集団の分散に近いことをサンプルサイズが 10 の場合について確かめるには次のようにする。なお，これまで何度か出てきたが，for ( ) { }という命令は，( ) 内の条件が満たされている間，{ }内の命令を繰り返せという意味の制御文である。

¶variance.R³

HT <- 144:191
NUM <- c(2,3,6,17,19,56,125,219,463,915,1609,2649,4550,7214,11005,16081,22098,29903,39048,48312,57703,66639,73332,78051,79829,77866,73767,66321,57993,48410,39081,29967,22055,15810,10875,7309,4596,2726,1519,939,462,224,128,50,31,14,5,4)
X <- rep(HT,NUM)
RNGkind("Mersenne-Twister",normal.kind="Inversion")set.seed(1)V <- rep(0,1000)
UV <- rep(0,1000)for (i in 1:1000) {s10 <- sample(X,10)
m10 <- mean(s10)V[i] <- sum((s10-m10)^2)/10
UV[i] <- sum((s10-m10)^2)/9}
layout(c(1,2))
hist(V,main="サイズ 10 の標本の分散の分布",xlim=c(0,80))
hist(UV,main="サイズ 10 の標本の不偏分散の分布",xlim=c(0,80))
print(c(mean(V),mean(UV)))µ´

この結果をみると，標本の分散も不偏分散も右裾を引いた分布になっていて，分布の位置がずれている。標本の分散の期待値よりも不偏分散の平均（つまり期待値）の方が母集団の分散に近いことが確認できる。

中心極限定理標本平均の分布については，もっと多くのことがいえる。実行するたびに標本抽出されるサンプルは異なるのだが，その平均は一定の法則に従うことが知られている。つまり，ほとんどの任意の母集団から抽出された√無作為標本の平均 X̄ の分布は，標本の大きさ n が増大するにつれて，平均 µ，標準偏差 σ/ n の正規分布に5近づく∗2 これを中心極限定理という。R で試してみよう。次のプログラムは 0 から 100 まで一様分布する乱数 5000 個と平均 100，標準偏差 10 の正規分布に従う乱数 5000 個をあわせた母集団（元の分布）から，標本の大きさ（サンプルサイズ）を増やしながら標本抽出するものである。

¶it04-2-2006.R³X <- c(runif(5000,0,100),rnorm(5000,100,10))
tsd <- function(X) { sqrt(var(X)*(length(X)-1)/length(X)) }layout(matrix(c(1,3,2,4),2,2))
hist(X,xlim=c(0,140),freq=F,main="元の分布")
Z5 <- rep(0,1000)for (i in 1:1000) { Z5[i] <- mean(sample(X,5)) }hist(Z5,xlim=c(0,140),freq=F)
Y2 <- dnorm(0:140,mean(X),tsd(X)/sqrt(5))
lines(0:140,Y2,col="red")
Z30 <- rep(0,1000)for (i in 1:1000) { Z30[i] <- mean(sample(X,30)) }hist(Z30,xlim=c(0,140),freq=F)
Y3 <- dnorm(0:140,mean(X),tsd(X)/sqrt(30))
lines(0:140,Y3,col="red")Z200 <- rep(0,1000)for (i in 1:1000) { Z200[i] <- mean(sample(X,200)) }hist(Z200,xlim=c(0,140),freq=F)
Y4 <- dnorm(0:140,mean(X),tsd(X)/sqrt(200))lines(0:140,Y4,col="red")µ´

サンプルサイズが大きくなるほど，正規分布に近づくと同時に，標本平均の標準偏差（つまり標準誤差）が小さくなっていき，ヒストグラムの幅が狭くなって，理論分布（赤い線）に近づくのが一目瞭然である。
前回やったように，生データのばらつきを示す指標として標準偏差が適切なのは，生データの分布が正規分布に近い場合に限られるが，中心極限定理から考えると，平均の信頼性を示す値としての標準誤差は，データの分布によらず，ある程度サンプルサイズが大きければ常に適切な指標となる。

信頼区間
標本平均の期待値は，母平均に一致するので，標本平均は，母平均のよい点推定量になっている。しかし，必ず一致するとは限らない。そこで，それがどれくらい確からしい推定かということを示すことを考える。通常，標本から計算される，ある区間の中に，正しく µ を含む割合が 95%であるような区間を推定したい（つまり，95%の信頼度をもって区間推定をしたい）と考える。
標本が大きい場合は，X̄ の正規分布において，ちょうど 95%の確率を囲むような最短の範囲を選択すればよい。そのために，正規分布の両側の裾 2.5%分を除く中央部分をとることができる。この部分は，母平均 µを中心として，大きい方と小さい方に，X̄ の標準偏差（つまり標準誤差）の 1.96 倍（1.96 は標準正規分布の97.5%点である）だけ動かした範囲を含む。つまり，
σσP r(µ − 1.96 √ < X̄ < µ + 1.96 √ ) = 0.95nn∗2 µ は母平均，σ は母集団の標準偏差である。

6である。この式を変形すると，

σσP r(X̄ − 1.96 √ < µ < X̄ + 1.96 √ ) = 0.95nnとなる。σ は未知なので，通常，標本の不偏標準偏差 sd(X) を用いる。つまり，µ の 95%信頼区間は，√√
X̄ − 1.96sd(X)/ n から X̄ + 1.96sd(X)/ n までとなる。

標本が小さい場合は，未知の σ に対して sd(X) を代入することが無視できない誤差の原因となる。そこで，

sd(X) を使って定義される t という統計量，t=X̄ − µ√sd(X)/ nが自由度 n − 1 の t 分布に従うことから，t 分布の 2.5%点から 97.5%点までを 95%信頼区間とすることで，この誤差を回避することができる。自由度が無限大になれば t 分布は正規分布に一致するので，結局，常にこちらで計算すればよいことになる。
すなわち，標本サイズ n，標本平均 X̄ ，標本の不偏標準偏差 sd(X) のとき，母平均の 95%信頼区間は，√X̄ − t0.025 sd(X)/ nから√X̄ + t0.025 sd(X)/ nまでになる。R で自由度 n − 1 の t 分布の 97.5%点を与える関数は qt(0.975,n-1) なので，R のプログラム
上で，例えば it04-1-2006.R で定義した X からのサイズ 100 のサンプル s100 から母平均の 95%信頼区間を推定するには，次のようにすればよい。

¶it04-3-2006.R³

HT <- 144:191
NUM <- c(2,3,6,17,19,56,125,219,463,915,1609,2649,4550,7214,11005,16081,22098,29903,39048,48312,57703,66639,73332,78051,79829,77866,73767,66321,57993,48410,39081,29967,22055,15810,10875,7309,4596,2726,1519,939,462,224,128,50,31,14,5,4)
X <- rep(HT,NUM)
RNGkind("Mersenne-Twister",normal.kind="Inversion")set.seed(1)s100 <- sample(X,100)
barX <- mean(s100)
sdX <- sd(s100)
t975 <- qt(0.975,length(s100)-1)
rootn <- sqrt(length(s100))print(barX - t975*sdX/rootn)
print(barX + t975*sdX/rootn)µ´

この結果から，95%信頼区間を表示するときは，[167.47, 169.19] のように記載するのが普通である。ただし，もちろん，R にはもっと楽にこの計算をしてくれる関数が用意されていて，標本データが付値されている
数値型変数 X について，t.test(X) とすれば，母平均がゼロであるという帰無仮説の検定結果（検定については後日説明する）とともに，母平均の推定値と母平均の 95%信頼区間が表示される。t.test(s100) が表示する母平均の推定値と 95%信頼区間は，上のコードで得られる結果と一致するはずである。各自確認されたい。

7

自由度
不偏分散の計算のところで出てきたが，自由度について直感的に理解するためには，次のように考えることもできる。n 個の数値からなる標本には，はじめ n の自由度がある。しかし，自由度の１つは，X̄ を計算するときに使われてしまい，sd(X) を計算するための偏差 Xi − X̄ に対しては n − 1 の自由度しか残らない。逆に見れば，サイズ n の標本について，はじめの n − 1 個の偏差は自由だが，最後の偏差は，全偏差の合計がゼロにならねばならない（そのように平均という量はとられている）ために決定済みであり，結局自由に決められる個数は n − 1 となる。
一般には，データの数から，推定した母数の数を引いた値が，その統計量や分布の自由度になる。

課題
50 から 99 までの値が 1000 ずつ，合計 50000 個の値からなる母集団があるとする。そこから 5 個のサンプルを 100 回抽出した時と 25 個のサンプルを 100 回抽出したときの標本平均の分布を図示してパワーポイントに貼り付け，2 つの分布を比較して考察したことを記入し，学生証番号と氏名を打って印刷したものに氏名を自筆して提出せよ。

8