医学情報処理演習第５回「データの分布と検定の概念」∗1
2005 年 11 月 7 日 中澤 港 (nminato@med.gunma-u.ac.jp)¶
前回の課題の回答例³(1) に回答した方はゼロだったので (1) は省略します。
(2) の回答例は，以下の通りです。
回答数は 311，標本平均が 1.385，標本の不偏分散は 0.0252 なので，母平均の 95%信頼区間の下限は，¶³

1.385 - qt(0.975,310)*sqrt(0.0252/310)µ´

より 1.367259 と得られ，上限は，¶³

1.385 + qt(0.975,310)*sqrt(0.0252/310)µ´

より 1.402741 と得られるので，95%信頼区間は [1.367, 1.403] である。
なお，小数以下何桁目まで表示するかという点に関しては，いくつかの考え方があるが，有効数字の１つ下の位で四捨五入してしまうか，正確な 95%信頼区間を間違いなく含む区間にするために，下限では切り捨て，上限では切り上げると考えるか，どちらかにすることが普通であるa 。
a この課題では，どちらの考え方でも回答は一致する。

µ´

はじめに
量的なデータの場合はヒストグラムや正規確率プロットによって分布の様子をみることができるし，カテゴリデータの場合は各カテゴリの度数分布図（あるいは割合）をみることによって，分布の様子をみることができる。今回は，いくつかの理論分布を紹介し，データの分布が理論分布に適合しているかどうかを調べる検定法を紹介し，併せて検定の考え方そのものにも説明を加える。仮説検定という考え方は，伝統的な統計解析の中では，かなり重要な部分を占めてきたので，ここできちんと整理しておく。

ベルヌーイ試行と２項分布
まずはカテゴリデータの分布から説明する。１回の実験で事象 S か事象 F のどちらかが起こり，しかもそれらが起こる可能性が，P r(S) = p, P r(F ) = 1 − p = q で何回実験しても変わらないとき，これを「ベルヌーイ試行」という。ベルヌーイ試行では，事象 F は事象 S の余事象になっている。
例えば，不透明な袋に黒い玉と白い玉が 500 個ずつ入っていて，そこから中を見ないで１つの玉を取り出して色を記録して（事象 S は「玉の色が黒」，事象 F は「玉の色が白」）袋に戻す実験はベルヌーイ試行である∗2 。
ベルヌーイ試行を n 回行って，S がちょうど k 回起こる確率は，P r(X = k) =n Ck pk q n−k である。n Ckは言うまでもなく n 個のものから k 個を取り出す組み合わせの数である。２項係数と呼ばれる。このような確率変数 X は，「２項分布に従う」といい，X ∼ B(n, p) と表す。E(X) = np, V (X) = npq である。
∗1 本資料は http://phi.med.gunma-u.ac.jp/medstat/it05.pdf としてダウンロード可能である。
∗2 注：袋に戻さないと１回実験するごとに事象の生起確率が変わっていくのでベルヌーイ試行にならない。なお，サンプリングとみれば，これは復元抽出である。

1

２項分布のシミュレーション
正二十面体（各面には１から 20 までの数字が割り振られている）サイコロを n 回 (n = 4, 10, 20, 50) 投げたときの，１から４までの目が出る回数を１試行と考えれば，これはベルヌーイ試行である。１回投げたときに１から４までの目が出る確率は理論的には 0.2（＝母比率は 0.2）と考えられるので，試行 1000 セットの度数分布を描く R のプログラムは下記の通り。

¶it05-1.R³times <- function(n) {dice <- as.integer(runif(n,1,21))
hit <- sum(ifelse(dice<5,1,0))
return(hit)}a <- c(4,10,20,50)
op <- par(mfrow=c(2,2))
for (i in 1:4) {y <- 1:1000
for (k in 1:1000) { y[k] <- times(a[i]) }barplot(table(y),main=paste("n=",a[i]))}par(op)µ´２項分布の理論分布
この例で，各 n についての理論的な確率分布は，P r(X = k) =n Ck 0.2k 0.8n−k なので，図を描くための Rのプログラムは下記の通り。

¶
it05-2.R³

op <- par(mfrow=c(2,2))
a <- c(4,10,20,50)
for (i in 1:4) {n <- a[i]k <- 0chk <- 1:(n+1)
names(chk) <- 0:n
while (k <= n) {chk[k+1] <- choose(n,k)*(0.2^k)*(0.8^(n-k))k <- k+1}
barplot(chk,main=paste("n=",n))}par(op)µ´

た だ し ，前 回 触 れ た dnorm() や dt() な ど ，R に は 様 々 な 確 率 分 布 に つ い て の 関 数 が あ り ，choose(n,k)*(0.2^k)*(0.8^(n-k)) は dbinom(k,n,0.2) と 同 値 で あ る 。こ の よ う に ，確 率 変 数 が取りうる各値に対して，その値をとる確率を与える関数を確率密度関数という。値が小さいほうからそれを全2部足した値を与える関数（つまり，その確率変数の標本空間の下限から各値までの確率密度関数の定積分）を分布関数（あるいは確率母関数，累積確率密度関数）と呼ぶ。
２項分布の確率変数の定義域は整数値なので，飛び飛びの値となる。その意味で，このような分布を離散分布という。離散分布には，２項分布の他には，ポアソン分布などがある。それに対して，正規分布や t 分布など，確率変数の定義域が実数である分布を，連続分布という。

正規分布
n が非常に大きい場合は，２項分布 B(n, p) の確率 P r(X = np + d) という値が，µ¶d21√exp −2npq2πnpqで近似できる。一般にこの極限（n を無限大に限りなく近づけた場合）である，µ¶(x − µ)2P r(X = x) = √exp −2σ 22πσ 21という形をもつ確率分布を正規分布と呼び，N (µ, σ 2 ) と書く。

z = (x − µ)/σ と置けば，1P r(Z = z) = √ exp2πµ−z 22¶となる。これを標準正規分布と呼び，N (0, 1) と書く。

R で標準正規分布の確率密度関数を [−5, 5] の範囲でプロットするには，curve(dnorm(x),-5,5) と打てば良い。curve() 関数は，連続分布をプロットするときに，定義域を x として，始点と終点をコンマで区切って与えれば曲線を描画してくれるので，とても便利な関数である。重ね描きする場合は，add=T を引数リストに加える。例えば，同じ範囲で平均 1，標準偏差 2 の正規分布の確率密度関数を赤い破線で重ね描きするには，

curve(dnorm(x,1,2),add=T,col="red",lty=2) とすればよい。
標準正規分布の 97.5%点（その点より小さい値をとる確率の積分値が 0.975 になるような点）を得るには，qnorm(0.975) とすればよいし，-1.96 より小さな値をとる確率を得るには，pnorm(-1.96) とすればよい。

χ2 分布
X1 , X2 , ..., Xv が互いに独立に標準正規分布 N (0, 1) に従うとき，
V =v
X

Xi 2i=12の分布を自由度 v の χ 分布という。この分布の確率密度関数は，

f (x|v) =³ v ´v/2−1³ v´
1exp −2Γ(v/2) 22である。
なお，言うまでもないが，Γ はガンマ関数で，正の実数 α に対して，Z ∞Γ(α) =xα−1 exp (−x)dx03であり，正の整数 α に対しては Γ(α) = (α − 1)! である。

E(x) = v で あ り ，V (x) = 2v で あ る 。自 由 度 1 の χ2 分布 を [0, 10] の範 囲で プ ロッ ト する に は，curve(dchisq(x,1),0,10) とすればよい。他の自由度のものを重ね描きするには，例えば自由度 2 のχ2 分布を赤破線で重ね描きしたければ，curve(dchisq(x,2),0,10,add=T,col="red",lty=2) とすればよい。
自由度 1 の χ2 分布の 95%点を得るには，qchisq(0.95,1) とすればよいし，3.84 より小さな値をとる確率を得るには，pchisq(3.84,1) とすればよい。

t 分布
標準正規分布に従う確率変数 U と，自由度 v の χ2 分布 χ2 (v) に従う確率変数 V があり，それらが独立のとき，

T = U/pV /vが従う分布のことをステューデントの t 分布という。この確率密度関数はf (t) = √µ¶−(v+1)/2
Γ((v + 1)/2)
t21+vvΓ(1/2)Γ(v/2)

である。これは，ステューデントというペンネームで論文を書いていたギネス社の技師ゴセット (Gosset WS)が初めて導いた分布である。
自由度 20 の t 分布の確率密度関数を [−5, 5] の範囲でプロットするには，curve(dt(x,20),-5,5) とすればよい。これが標準正規分布より裾が長い分布であることを見るために標準正規分布を赤い点線で重ね描きするには，続けて，curve(dnorm(x),-5,5,add=T,col="red",lty=2) とすればよい。
また，自由度 20 の t 分布の 97.5%点を得るには，qt(0.975,20) とすればよいし，2 より小さな値をとる確率を得るには，pt(2,20) とすればよい。

検定の考え方と第一種，第二種の過誤検定とは，帰無仮説（一般には，差がない，という仮説）の元で得られた統計量を，既知の確率分布をもつ量と見た場合に，その値よりも外れた値が得られる確率（これを「有意確率」と呼ぶ）がどれほど小さいかを調べ，有意水準∗3 より小さければ，統計的に意味があることと捉え（統計的に有意である，という），帰無仮説がおかしいと判断して棄却する（つまり，「差がないとは言えない」と判断する）という意思決定を行うものである。
この意思決定が間違っていて，本当は帰無仮説が正しいのに，間違って帰無仮説を棄却してしまう確率は，有意水準と等しいので，その意味で，有意水準を第一種の過誤と（αエラーとも）呼ぶ（逆に，本当は帰無仮説が正しくないのに，その差を検出できず，有意でないと判断してしまう確率を，第二種の過誤と（βエラーとも）呼び，1 から第二種の過誤を引いた値が検出力になる）。
注意しなければいけないのは，有意確率の小ささは，あくまで，帰無仮説のありえなさを示すだけであっ∗3 分析者が決める一定の確率。当該研究分野の伝統に従うのが普通である。先行研究があればそれに従う。他に基準がなければ5%か 1%にすることが多い。

4

て，差の大きさを意味するのではない点である。ここを勘違いしたレポートなどが時折見られるので注意されたい。

両側検定と片側検定
２つの量的変数 X と Y の平均値の差の検定をする場合（平均値の差の検定についてはまた次回詳しく触れる），それぞれの母平均を µX ，µY と書けば，その推定量は µX = mean(X) =PPX/n と µY = mean(Y ) =Y /n となる。

両側検定では，帰無仮説 H0 : µX = µY に対して対立仮説（帰無仮説が棄却された場合に採択される仮説）

H1 : µX 6= µY である。H1 を書き直すと，「µX > µY または µX < µY 」ということである。つまり，t0 を「平均値の差を標準誤差で割った値」として求めると，t0 が負になる場合も正になる場合もあるので，有意水準 5%で検定して有意になる場合というのは，t0 が負で t 分布の下側 2.5%点より小さい場合と，t0 が正で t 分
布の上側 2.5%点（つまり 97.5%点）より大きい場合の両方を含む。t 分布は原点について対称なので，結局両側検定の場合は，上述のように差の絶対値を分子にして，t0 の t 分布の上側確率∗4 を２倍すれば有意確率が得られることになる。
片側検定は，先験的に X と Y の間に大小関係が仮定できる場合に行い，例えば，X の方が Y より小さくなっているかどうかを検定したい場合なら，帰無仮説 H0 : µX ≥ µY に対して対立仮説 H1 : µX < µY となる。この場合は，t0 が正になる場合だけ考えればよい。有意水準 5%で検定して有意になるのは，t0 が t 分
布の上側 5%点（つまり 95%点）より大きい場合である。なお，R で平均値の差の検定を行うための関数はt.test() だが，詳しくは次回説明する。

分布の正規性の検定
高度な統計解析をするときには，データが正規分布する母集団からのサンプルであるという仮定を置くことが多いが，それを実際に確認することは難しいので，一般には，分布の正規性の検定を行うことが多い。考案者の名前から Shapiro-Wilk の検定と呼ばれるものが代表的である。

Shapiro-Wilk の検定の原理をざっと説明すると，Zi = (Xi − µ)/σ とおけば，Zi が帰無仮説「X が正規分布にしたがう」の下で N (0, 1) からの標本の順序統計量となり，c(i) = E[Z(i)], dij = Cov(Z(i), Z(j)) が母数に無関係な定数となるので，「X(1) < X(2) < ... < X(n) の c(1), c(2), ...c(n) への回帰が線形である」を帰無仮説として，そのモデルの下で σ の最良線形不偏推定量 ˆ(σ) =PnPn22
i=1 (Xi − X̄) をi=1 ai X(i) と S =Pn用いて，W = (kσ̂ 2 )/S 2 を検定統計量として検定するものである。k は i=1 (kai )2 = 1 より求められる。

R で数値型変数 X の分布が正規分布にフィットしているかどうかを検定するには，¶³shapiro.test(X)µ´とすればよい。変数 X のデータ数（ベクトルの要素数，R のコードでは length(X)）は，3 から 5000 の間でなければならない。2 以下では分布を考える意味がなく，また，検定統計量 W の分布がモンテカルロシミュレーションによって得られたものであるため，あまりに大きなサンプルサイズについては値が与えられて∗4 t 分布の確率密度関数を t0 から無限大まで積分した値，即ち，t 分布の分布関数の t0 のところの値を 1 から引いた値。R では

1-pt(t0, 自由度)。

5

いないのである。
正規性の検定にはたくさんの方法が提案されているが，ここではもう１つだけ紹介しておこう。提案者の名前から，ギアリー (Geary) の検定と呼ばれるものである。現在のところ，R にはデフォルトでは入っていな
いが，比較的簡便で使いやすい検定である。以下，簡単に原理を説明する。
左右対称な分布について，裾の長さを，平均値のまわりの１次の絶対モーメントを２次のモーメントの平方根で割ったもので測ることにすると，その一致推定量 G が，Pni=1 |Xi − X̄|
G= q Pnn i=1 (Xi − X̄)2となる。この G を用いて帰無仮説 H0:「データ X が正規分布からの標本」を検定することができる。対立仮説の下での分布が正規分布よりも裾が長い対称分布（例えば t 分布のような）であれば G < g0 のとき帰無仮説を棄却する。g のパーセント点については，uα を標準正規分布の 100α%点として，n が大きければ近似的にrg(α; n) '2+ uαπr1−3 1√π nで得られることがわかっている。そこで，R のプログラムとしては下枠内のようにして，例えば G がg(0.05) の値よりも小さければ，正規分布よりも 5%水準で有意に裾が長いとして帰無仮説を棄却すればよい。下枠内の関数定義をした後で，試しに geary.test(rnorm(1000)) とすれば，ほぼ全ての試行で G はg(0.05) よりも大きくなるし，geary.test(rt(1000,20)) とすれば，ほぼ全ての試行で G は g(0.01) よりも小さくなるであろう。

¶it05-3.R³

geary.test <- function(X) {m.X <- mean(X)l.X <- length(X)G <- sum(abs(X-m.X))/sqrt(l.X*sum((X-m.X)^2))
g5 <- sqrt(2/pi)+qnorm(0.05)*sqrt(1-3/pi)/sqrt(l.X)
g1 <- sqrt(2/pi)+qnorm(0.01)*sqrt(1-3/pi)/sqrt(l.X)
cat("G=",G," / g(0.05)=",g5," / g(0.01)=",g1,"\n")}µ´

なお，作図の説明で触れた qqnorm(X) をしてみるのも，分布の正規性をチェックするにはいい方法である。
正規分布に従っていれば直線に乗るはずであり，外れているときにどのように外れているかが見える。

課題
次の (1)(2) のいずれか（両方でもよい）に答えてください。結果は配布する紙に学籍番号，氏名と共に自筆して提出すること。結果の提出をもって出席確認とします。

(1)

http://phi.med.gunma-u.ac.jp/medstat/p01.txt を読み込み，男女（変数 SEX で与えられる）別に，体重（変数 WT で与えられる）が正規分布に従っていると言えるかどうか検定してください。答えだけでなく手順も書いてください。

(2)

今回も含めて，これまでの演習への意見・質問・感想を書いてください。

6