sec6

８．仮説検定（Test of hypothesis）

帰無仮説（H₀）と対立仮説（H₁）

統計学で扱う仮説（hypothesis）とは，母集団に対する断定や推測．たとえば，

母集団は正規分布に従っている．
母集団平均は 0 である．
母集団 A と母集団 B の平均は等しい．

などである．

統計的仮説検定で用いられる仮説は，まず，帰無仮説（null hypothesis）という形式で与えられる．
帰無仮説は棄却（reject）されることに意味がある仮説である．
帰無仮説と反対の仮説を対立仮説（altanative hypothesis）という．

上の3番目の例でみると，

帰無仮説：母集団 A と母集団 B の平均は等しい．（H₀： μ_A ＝ μ_B）

対立仮説：母集団 A と母集団 B の平均は等しくない．（H₁： μ_A ≠ μ_B）

母集団 A と母集団 B は異なる処理（薬の投与など）をしているので，実験の目的は，母集団 A と母集団 B の平均は異なる（処理効果がある）ことを言いたい（対立仮説が正しいことを望む）のだが，まずは「等しい（処理効果無し）」と仮定してみようという考え方で，数学の背理法と似た論理である．

背理法：√2 が無理数であることを証明するため，まず√2 が有理数であると仮定し，矛盾があることを示す．つまり，有理数であることは絶対ありえない（確率 0 である！）ことを示す．この矛盾は，そもそも√2 を有理数とした仮定が誤っていたからであると考え，有理数という仮定を棄却して，無理数であることを証明する．

検定の概要

検定統計量（Test statistic）

　標本から算出される量で，検定に用いられるもので，t 値（t value），F 値（F value）などがある．この値から帰無仮説を受託（accept）（採択）するか 棄却（reject）（対立仮説の採択）するかを判定する．

p 値（p value）と有意水準（Significance level）

　統計的仮説検定では，たとえば２つの母集団平均が等しいという帰無仮説を考えると，この帰無仮説のもとで，検定統計量（標本平均の差に基づく t 値など）以上（もしくは未満）の値が得られる確率を求める．R ではこの値が p 値で表示される．
　p 値はくだけた言い方をすれば，帰無仮説が正しいとしたときに，標本のようなデータが得られる確率である．
　これが十分小さい（ほとんどありえない）ときは，平均が等しいと仮定したことが誤りであったと判断して帰無仮説を棄却し，２つの母集団平均には差があると結論づける．
　この確率がそれほど小さくない場合は，このような統計量が得られることもありえると考え，帰無仮説を採択し，平均が等しいと考えてもよいとする．
　棄却か採択かの判断の基準となる確率を有意水準といい， 5 % や 1 % がよく用いられる．

　統計のソフトが発達していなかった頃は，検定統計量である t 値や F 値を電卓等で算出し，その値を t 分布や F 分布の 5 ％や 1 ％の有意水準に対応する数表と照らし合わせて検定を行い，5 ％有意とかを記述していた．
　現在では，ソフトが検定統計量に対する p 値を直接計算してくれるので数表はいらなくなった．この結果，検定に重要な数字であった t 値や F 値より，より直接的な p 値が重要な指標になってきた．p 値をみれば，何％有意かが一目でわかるので，結果にわざわざ 5 ％有意とかを記述する必要がなくなってきており，論文の書き方も変わってきている．

片側検定と両側検定

実験状況によっては，薬投与などの処理を行った集団（処理群）平均 μ_A が，薬を投与しない集団（対照群）の平均 μ_B より小さくなることはないことが事前にわかっているような場合がある．このようなとき，

帰無仮説，H₀： μ_A ＝ μ_B

対立仮説，H₁： μ_A ＞ μ_B

となる．これは，事前情報より，μ_A ＜ μ_B となる可能性をまったく考えない場合である．
このため検定には，片側 5 ％点や 1 ％点を用いる．

両側検定と信頼区間

母集団平均に対する両側検定は，母集団平均に対する信頼区間と大きな関係がある．いま，帰無仮説（H₀）と対立仮説（H₁）が，

H₀：μ ＝ μ₀
H₁：μ ≠ μ₀

であり，母分散 σ² が既知のときを考える．
標本の大きさがnで，標本平均が $\bar{x}$ であったとすると，母平均μに対する 95％信頼区間は，

${\rm Pr} \Bigl[ -1.96 < \frac{\sqrt{n}(\bar{x}-\mu)}{\sigma} < 1.96 \Bigr] = 0.95$

${\rm Pr} \bigl[ \bar{x} - 1.96 \cdot\frac{\sigma}{\sqrt{n}} < \mu < \bar{x} + 1.96 \cdot\frac{\sigma}{\sqrt{n}} \Bigr]$

となる．
　一方，この検定の検定統計量は，標本平均の標準化値の絶対値

$|z| =\frac{\sqrt{n} |\bar{x} - \mu_0 |}{\sigma}$

で，有意水準 5 ％で帰無仮説を受諾するのは，検定統計量 |z| が両側 5 ％点である 1.96 以下のときである．つまり，

帰無仮説を受諾 ⇔

$-1.96 < \frac{\sqrt{n}(\bar{x}-\mu)}{\sigma}} <1.96$

である．この両者の関係より，
帰無仮説を受諾 ⇔ 母平均の信頼区間に μ₀ が含まれる．

帰無仮説を棄却 ⇔ 母平均の信頼区間に μ₀ が含まれない．
が成り立つ

検定における２種類の過誤

検定は，仮説を棄却するか採択するかのいずれかであるが，統計量は分布をもつので，この判定には間違いが起こることがある．
以下のように，この過誤には 2 種類がある．

統計的検定における2種類の過誤
仮説の棄却（reject）仮説の採択（accept）

仮説が真（true）のとき第１種の過誤（Type 1 error）正解

仮説が偽（false）のとき正解第２種の過誤（Type 2 error）

統計的検定における2種類の過誤
	仮説の棄却（reject）	仮説の採択（accept）
仮説が真（true）のとき	第１種の過誤（Type 1 error）	正解
仮説が偽（false）のとき	正解	第２種の過誤（Type 2 error）

第１種の過誤が有意水準である．また，第２種の過誤の確率を β としたとき，仮説が偽のとき正しく仮説を棄却する確率，1 - β，を検出力，もしくは検定力（power）という．よい検定は，第１種の過誤を固定したもとで検出力の高い検定方式である．

正規母集団の母平均に対する t 検定

１つの母集団に対する検定（One sample problem）

　平均 μ，分散 σ² がともに未知である正規母集団に対して，

帰無仮説 H₀： μ ＝ μ₀

対立仮説 H₁： μ ≠ μ₀

の両側検定を考える．　いま，母集団から大きさ n の無作為標本 x₁，x₂，…，x_n を抽出したところ，標本平均が x^-，標本分散が s² であったとする．帰無仮説のもとでは，標本平均は，

と分布するので，これを，標本平均の標準誤差 s/√n で標準化した t は，

のように自由度 n - 1 の t 分布に従う．この分布の97.5％分位点を t(n - 1)_0.975 とすると，有意水準 5 ％の検定は，

|t| > t(n - 1)_0.975

のとき帰無仮説を棄却する．|t| が検定統計量で，この値を |t| 値という．

なお，この検定は，対のある標本に適用できる．対のある標本とは，n 組のペアー標本（paired smple），

(x₁，y₁ )，(x₂，y₂ )， …，(x_n，y_n )

からなる．正規性の仮定のもとでは，

x_i ～ N( μ_i，σ_x² )， y_i ～ N( μ_i ＋ δ，σ_y² )

ここで興味ある母数は δ であり，μ₁，…，μ_n は攪乱母数（nuisance parameter）である．y_i と x_i の差を取ると，

z_i ＝ y_i － x_i → z_i ～ N( δ，σ_z² )

となるので，１つの母集団に対する検定に帰着する．なおこの問題は，反復のない 2×n の２元配置と考えて解くこともできる．

園芸療法前後での MSE 得点
　園芸療法（horticultural therapy）は第2次大戦後，兵士のリハビリテーション等を通して欧米で発達してきた．日本でも90年代から身体障害者，精神障害者，高齢者等の生活の質（QOL, Quality Of Life）やセルフエスティーム（self-esteem，自尊感情）を向上させる手段として注目されてきていている．
　園芸療法の数値的･客観的な効果の計測を行うため，認知症診断に用いられる MSE（Mental State Examination，心理機能検査），CDT（Clock Drawing Test，時計描画テスト）と老年者や認知症患者の日常生活の遂行能力を測る NM スケール（N 式老年者精神状態評価尺度）， N-ADL（N 式老年者日常生活動作能力評価尺度），および，気分の状態（抑うつ，活気のなさ，怒り，疲労，緊張，混乱）を計測する POMS（Profile Of Mood State），自己記入式 QOL 質問表 QUIK を用いた．これらの調査を園芸療法の前後でクライアントに対して行い，効果の数値的計測を行った．
　ここでは，MSE の結果についての解析を行う．MSE は 30 点満点で，点数が高いほど日常生活に適応していると考えられている．MSE に関しては前後とも回答したクライアントの総数は 19 例であった．データは下の表にまとめられている．このデータを用いて，園芸療法を施す前と後では，心理機能に差がないという帰無仮説，すなわち，

H₀：差＝ 0

の検定を行う．

データダウンロード

文献

鈴木修（編），大森宏，児玉良治，渡辺俊之，矢野広，山根健治（2004）．専修学校における園芸療法士教育育成システムの研究開発（文部科学省委託平成15年度「専修学校先進的教育研究開発事業」）
鈴木修，渡辺俊之，矢野広，山根健治，大森宏，伊東正信，最上正秀，山下容子，小泉力，児玉良治，頭士智美，細井薫，水口聡子，遠藤久子，樋田奈穂子，小島ユリ，郡司敏幸（2005）．福祉サービス提供者に対する園芸療法教育システムの研究開発（平成16年度文部科学省「専修学校社会人キャリアアップ教育推進事業」）

# １つの母集団に対する t 検定の R スクリプト
engei <- read.csv("engei.csv") # csv データ読み込み　
engei # データの表示　
x <- engei[,3] # 園芸療法前の MSE 得点を x に格納（列数で指定）　
x <- engei$療法前 # 園芸療法前の MSE 得点を x に格納（変数名で指定）　
y <- engei$療法後 # 園芸療法後の MSE 得点を x に格納（変数名で指定）　
t.test(y, x, paired=TRUE) # １母集団 t 検定，paired = TRUE で対標本を指定
d <- y - x # 療法後－療法前で療法の効果をみる　
d # 療法の効果の表示　
t.test(d) # １母集団 t 検定，先ほどと同じ検定　
n <- length(d) # 標本の大きさ（サンプルサイズ）　
mean(d) # 標本平均　
sd(d) # 標本標準偏差　
dv <- n - 1 # 標本の自由度　
t <- sqrt(n)*mean(d)/sd(d) # 効果がないとの帰無仮説のもとでの t 値　
t # 検定統計量 t 値の表示　
2*(1 - pt(t, df=dv)) # 両側検定の p 値　
t0 <- qt(0.975, df=dv) # 両側 5 ％検定の閾値　
dw <- t0*sd(d)/sqrt(n) # 95％信頼区間の幅　
mean(d)-dw # 95％信頼区間の下限　
mean(d)+dw # 95％信頼区間の上限　
t.test(d, alternative="greater") # 片側検定　

　このデータでは，クライアント全体では園芸療法の効果は認められなかった．また，MSE の満点が30点なので，元々心理機能に問題がなく満点近い得点のクライアントでは，効果が認められないのも当然といえる．その上，園芸療法処方の前後で半年程度のタイムラグがあるので，何もしなくても認知症の症状が進行する場合もあり，効果が検出しずらいこともある程度予測できる．

課題：クライアントのうち若年層（ 85 歳以下）を取り出し，園芸療法の効果を検定せよ．

t 検定の有意水準
　t 検定の有意水準をシミュレーションにより確認してみる．平均 μ = 10，分散 σ² = 4 の正規分布から大きさ n = 10 の標本を抽出し，その標本平均 x^- と標本（不偏）分散 s² から t 値を計算する．これを N = 10000回繰り返すと，t 値のヒストグラムができる．このヒストグラムは自由度 9 の t 分布に従っているはずである．

帰無仮説，H₀：μ = 10

の検定では， t 値の絶対値 |t| が自由度 9 の t 分布の97.5％分位点 t(9)_0.975 より大きくなった場合に棄却される．このため，この個数を数えれば第１種の過誤の確率，つまり有意水準がシミュレーションにより推定できる． N = 10000回のシミュレーションでは，これがほぼ 5 ％になっていた．

**#t 検定の有意水準の R スクリプト**
N <- 10000	# シミュレーション回数
n <- 10	# サンプルサイズ
n1 <- rnorm(N*n, mean=10, sd=2)	# N(10, 4) から大きさ n の標本を N 回シミュレーション
n1.mat <- matrix(data=n1, ncol=n)	# データ行列
n1.mean <- apply(n1.mat, 1, mean)	# 各サンプルの平均
n1.var <- apply(n1.mat, 1, var)	# 各サンプルの分散
n1.var3 <- n1.var/n	# 各標本平均の分散
n1.td <- (n1.mean - 10)/sqrt(n1.var3)	# 各サンプルの t 値
mean(n1.td)	# t 値の平均
sd(n1.td)	# t 値の標準偏差
mt <- ceiling(max(abs(n1.td)))	# t 値の絶対値の最大値
xq1 <- qt(0.025, df = (n-1))	# t(9) の 2.5% 点
xq2 <- qt(0.975, df = (n-1))	# t(9) の 97.5% 点
length( n1.td[abs(n1.td) > xq2])	# 5% 検定で有意となった個数
hist(n1.td, breaks=seq(-mt, mt, by=0.2), probability=TRUE, xlab="t 値", ylab="密度", main="")
abline(v=xq1, col="red")	# 採択域の下限
abline(v=xq2, col="red")	# 採択域の上限
#curve(dt(x, df=(n-1)), -6, 6, add=T, col="red")	# t(9) の表示
title(main="t 値のヒストグラムと有意となった個数")	# タイトル

非心 t 分布（non-central t distribution）
　帰無仮説が正しいときに，t 値の従う分布は t 分布に従うことを見た．帰無仮説が正しくないときに t 値の従う分布を非心 t 分布という．いま，母集団分布を N(μ，σ²) としたとき，

帰無仮説：H₀：μ = μ₀，対立仮説：H₁：μ ≠ μ₀

の t 検定を，大きさ n の標本で行う場合を考える．このとき標本平均が x^-，標本標準偏差が s であったとすると，検定統計量

t = √n(x^- - μ₀)/s

の分布は，帰無仮説が正しくないとき，自由度 n - 1，非心パラメーター（non-central parameter）

δ = √n(μ - μ₀)/σ

の非心 t 分布に従う．この分布を用いると，t 検定の検出力を求めることができる．

**# 非心 t 分布の R スクリプト**
n <- 10 # サンプルサイズ m0 <- 10 # 帰無仮説の平均 s <- 2 # 標準偏差　 m <- 11 # 真の平均 delta <- sqrt(n)*(m-m0)/s # 非心度 curve(dt(x, df=(n-1)), -4,7) curve(dt(x, df=(n-1), ncp=delta), -4,7, col="red", add=T) title(main="non-central t distribution")

# 非心 t 分布の R スクリプト

n <- 10 		# サンプルサイズ 
m0 <- 10		# 帰無仮説の平均
s <- 2		# 標準偏差　
m <- 11		# 真の平均 
delta <- sqrt(n)*(m-m0)/s	# 非心度
curve(dt(x, df=(n-1)), -4,7)
curve(dt(x, df=(n-1), ncp=delta), -4,7, col="red", add=T)
title(main="non-central t distribution")

t 検定の検出力（power）
　前節では帰無仮説が真のとき，t 検定により仮説が棄却する確率を求めた．では，帰無仮説が偽であるとき，帰無仮説がどれくらいの確率で正しく棄却できるであろうか．すなわち，t 検定の検出力はどの程度であろうか．この問題に解析的に答えるのは難しいが，シミュレーションなら簡単に求めることができる．
　いま，母集団の平均を m = 11 として，N(11，4) から大きさ n = 10 の標本を抽出する．帰無仮説が偽であるとき，標本から得られる t 値の分布は簡単ではないが，N = 10000回のシミュレーションを行えば t 値のヒストグラムが得られる．シミュレーションで得られた t 値の中で帰無仮説が棄却される割合を求めれば検出力が推定できる．
　シミュレーションの結果では，検出力は約 30 ％であった．

**#t 検定の有意水準の R スクリプト**
m <- 11	# 母集団平均（帰無仮説が偽の場合）
n1 <- rnorm(N*n, mean=m, sd=2)	# N(m, 4) から大きさ n の標本を N 回シミュレーション

課題１：：母集団平均を m = 12 としたときの検出力を求めよ．

課題２：：真の母集団平均が μ = 11 であるとき，検出力を 90％程度にするために必要なサンプルサイズを求めよ．

２つの母集団に対する検定（Two sample problem）

　２つの母集団 A，B があり，それぞれが平均を μ_A，μ_B，分散を σ_A²，σ_B² の正規分布に従っているが，その値は未知であるとする．いま，両集団の分散の値が等しく， σ_A²＝σ_B²＝σ²，と仮定できるとしよう．このとき，２つの母平均に対する

帰無仮説，H₀： μ_A ＝ μ_B

対立仮説，H₁： μ_A ≠ μ_B

の検定は t 分布を用いて行える．

　母集団 A から大きさ n_A，母集団 B から大きさ n_B の標本を抽出した．母集団 A からの標本の標本平均が x^-_A，標本分散が s_A² であり，母集団 B の標本平均が x^-_B，標本分散が s_B² であった．母集団 A，B が共通の分散 σ² をもつとすると，その推定値 s² は以下のように推定できる．

母集団 A からの標本の偏差平方和：　S_A＝（n_A－1）s_A²

母集団 B からの標本の偏差平方和：　S_B＝（n_B－1）s_B²

母集団 A，B 全体での偏差平方和：　S ＝ S_A ＋ S_B ＝（n_A－1）s_A²＋（n_B－1）s_B²

母集団 A，B 共通の標本分散： bunsan

また，母集団Aの標本分布は，N（μ_A，σ²）であり，母集団Bでは， N（μ_B，σ²）であることから，それぞれの標本平均は，

x^-_A ～ N（μ_A，σ²/n_A）， x^-_B ～ N（μ_B，σ²/n_B）

と分布する．これより，標本平均の差x^-_A－ x^-_Bは，

と分布する．

帰無仮説（H₀： μ_A ＝ μ_B）のもとでは，μ_A－μ_B＝0，なので，標本平均の差は，

と分布する．これを標準化した z 値，

において，標準偏差 σ の代わりに標本標準偏差 s を代入した t 値が，

のように自由度 n_A＋n_B－2 の t 分布に従うことを利用して検定ができる．なお，母集団 A，B からの標本の大きさがともに等しく， n_A＝n_B＝n であるときは，式がずっと簡単になる．

母集団A，Bで分散の同等性が疑われるときは，ウェルチ（Welch）の検定を用いる．

園芸療法区と対照区の MSE 得点
　前節では，心理機能検査（MSE）データを対のある 2 標本データとして解析した．この節では，このデータの園芸療法前データを対照とし，園芸療法後データを処理とみなして解析してみる．この場合，実際とは異なりこのデータの解析法としては正しくないが，対照 19 名，処理 19 名の計 38 名のデータとみなした場合である．　まず最初に，処理区と対照区での分布の様子を，箱ひげ図（box plot）で見て両者の違いの程度を観察し，その後２つの母集団に対する t 検定（２標本 t 検定）を行う．

**# ２標本 t 検定の R スクリプト**
engei <- read.csv("engei.csv")	# csv データ読み込み
x <- engei$療法前	# 園芸療法前の MSE 得点を x に格納（変数名で指定）
y <- engei$療法後	# 園芸療法後の MSE 得点を y に格納（変数名で指定）
boxplot(x, y, names= c("対照区","処理区"), ylab="MSE 得点", cex.axis=0.8)
title(main="園芸療法処理区と対照区の MSE 得点分布の箱ひげ図", cex.main=1.0)
t.test(y, x, var.equal=T)	# ２標本 t 検定
mx <- mean(x)	# 対照区平均
my <- mean(y)	# 処理区平均
nx <- length(x)	# 対照区サンプルサイズ
ny <- length(y)	# 処理区サンプルサイズ
dfv <- nx+ny-2	# データの自由度
vx <- var(x)	# 対照区標本分散
vy <- var(y)	# 処理区標本分散
v <- ((nx-1)vx + (ny-1)vy)/dfv	# こみにした共通の分散
d <- my - mx	# 処理区平均－対照区平均
vd <- v/nx + v/ny	# 処理平均の差の分散
tv <- d/sqrt(vd)	# t 値
2*(1 - pt(abs(tv), df=dfv))	# p 値
t0 <- qt(0.975, df=dfv)	# 自由度 dfv の t 分布の97.5％点
dw <- t0*sqrt(vd)	# 95 ％信頼区間の幅
d - dw	# 平均の差の 95％信頼区間の下限
d + dw	# 平均の差の 95％信頼区間の上限
t.test(x, y)	# ウェルチ（Welch）の検定

課題：クライアントのうち若年層（ 85 歳以下）を取り出し，２標本 t 検定を行え．

正規母集団の母分散に対する検定

母集団分散の検定

　　平均 μ，分散 σ² がともに未知である正規母集団に対して，

帰無仮説 H₀： σ² ＝ σ₀²

対立仮説 H₁： σ² ≠ σ₀²

の検定を考える．　いま，母集団から大きさ n の無作為標本 x₁，x₂，…，x_n を抽出したところ，標本平均が x^-，標本分散が s² であったとする．すると，帰無仮説ももとで（under H₀），標本分散に関係した量が，

と自由度 n - 1 の χ² 分布に従うので，U を検定統計量にして検定が行える．

　有意水準 5 ％の検定は，自由度 n - 1 の χ² 分布の 2.5％点と 97.5％点をそれぞれ χ²(n - 1)_0.025， χ²(n - 1)_0.975 とすると，

U ＜ χ²(n - 1)_0.025， U ＞ χ²(n - 1)_0.975

のいずれかの不等式を満たしたとき帰無仮説を棄却し，母分散は σ₀² と有意に異なると結論づける．

２つの母集団分散の同等性の検定

　２つの母集団 A，B があり，それぞれが平均を μ_A，μ_B，分散を σ_A²，σ_B² の正規分布に従っているが，その値は未知であるとする．このとき，２つの母分散の同等性の検定，

帰無仮説，H₀： σ_A² ＝ σ_B²

対立仮説，H₁： σ_A² ≠ σ_B²

の検定を考える．

　母集団 A から大きさ n_A，母集団 B から大きさ n_B の標本を抽出した．母集団 A からの標本の標本平均が x^-_A，標本分散が s_A² であり，母集団 B の標本平均が x^-_B，標本分散が s_B² であるとする．すると，標本分散に関係した量はそれぞれ

と χ² 分布に従い，それぞれが独立である．これらの量の比は，

のように，自由度 n_A - 1，n_B - 1 の F 分布に従う．

　ところで，帰無仮説が正しいとすると，σ_A² ＝ σ_B² とおけるので，母集団の分散比は，γ₀ = σ_A²/σ_B² = 1，となる．このとき，標本分散の分散比の統計量 γ が，

と，自由度 n_A - 1，n_B - 1 の F 分布に従うので，この γ 値を検定統計量にして２つの母分散が等しいという帰無仮説の検定が行える．
　すなわち，有意水準 5 ％の検定を行うには，自由度 n_A - 1，n_B - 1 の F 分布の 2.5％点と 97.5％点をそれぞれ F(n_A - 1，n_B - 1)_0.025， F(n_A - 1，n_B - 1)_0.975 とすると，検定統計量 γ が，

γ ＜ F(n_A - 1，n_B - 1)_0.025， γ ＞ F(n_A - 1，n_B - 1)_0.975，

のいずれかの不等式を満たしたとき帰無仮説を棄却し，２つの母集団の分散は有意に異なると結論づける．

２つの母集団の分散比の信頼区間

　２つの母集団 A，B の分散 σ_A²，σ_B² の分散比， γ₀ = σ_A²/σ_B²，の 95％信頼区間は，上記の考えから簡単に求めることができる．すなわち，互いに独立に χ² 分布する変量の比が，標本分散の分散比 γ と母集団分散比 γ₀ の比となり，

と分布する．これより，母集団分散比の 95％信頼区間は，

となる．

園芸療法区と対照区の MSE 得点
　前節では，心理機能検査（MSE）データを，検定の勉強に用いるため，対照区と処理区とみなして平均の比較をおこなった．両区の分散の等質性が疑われるので，分散の等質性の検定を行ってみる．

# 分散等質性の R スクリプト
engei <- read.csv("engei.csv") # csv データ読み込み　
x <- engei$療法前 # 園芸療法前の MSE 得点を x に格納（変数名で指定）　
y <- engei$療法後 # 園芸療法後の MSE 得点を y に格納（変数名で指定）　
var(x) # 対照区の標本分散　
var(y) # 処理区の標本分散　
var.test(x, y) # 分散の等質性の検定　

課題：分散の等質性検定で表示される数値を R の基本統計関数を用いて出せ．

分散比 F 検定の検出力
　２つの母集団の分散の違いはどれくらいあれば検出可能であるか，シミュレーションで調べる．
　まず，真の分散比が１のときを N = 10000 回シミュレーションし，正しい有意水準が得られるかを確認する．次に，真の分散比が２のときの検出力を２つの母集団からそれぞれ n = 20 のサンプルを抽出して行う．
分散の推定精度は標本の大きさ（サンプルサイズ）の影響を大きく受けると考えられるので，いくつかのサンプルサイズで試す．

**#分散比 F 検定の検出力の R スクリプト**
N <- 10000	# シミュレーション回数
n <- 20	# サンプルサイズ
s1 <- 1; s2 <- 2	# 母集団分散（分散比 = 2）
n1 <- rnorm(N*n, mean=0, sd=sqrt(s1))	# 正規乱数
n2 <- rnorm(N*n, mean=0, sd=sqrt(s2))	# 正規乱数
n1.mat <- matrix(data=n1, ncol=n)	# N×n データ行列（母集団１）
n1.var <- apply(n1.mat, 1, var)	# 各行の分散
n2.mat <- matrix(data=n2, ncol=n)	# N×n データ行列（母集団２）
n2.var <- apply(n2.mat, 1, var)	# 各行の分散
ratio <- n2.var/n1.var	# 分散比
f0 <- qf(0.025, df1=(n-1), df2=(n-1))	# F 分布2.5％点
f1 <- qf(0.975, df1=(n-1), df2=(n-1))	# F 分布97.5％点
m0 <- length(ratio[ratio < f0])	# 2.5％以下の個数
m1 <- length(ratio[ratio > f1])	# 97.5％以下の個数
m0+m1	# 帰無仮説を棄却した個数
mf <- ceiling(max(ratio))	# 分散比の最大値
hist(ratio, breaks=seq(0, mf, by=0.2), probability=TRUE, xlab="分散比", ylab="密度", main="")
abline(v=f0, col="red")	# F 分布2.5％点の表示
abline(v=f1, col="red")	# F 分布97.5％点の表示
title(main="分散比のヒストグラムと有意となった個数")	#

課題：淡水性ウナギの汽水域での生理活性の違いのデータ、ウナギデータ，を２標本データとして違いを検定せよ．ただし，２母集団で分散の違いがあるようなので，まず，分散の同等性の検定を行い，違いが認められるときは，データの対数変換を行うなどして，分散安定化変換を行って t 検定を行ってみよ．また，データの対数変換などを行わず．ウェルチ検定やノンパラメトリック検定を行った場合との比較検討を行え．

成功確率（比率）に関する検定

標準正規分布による近似検定（大標本理論）

　成功確率 p のベルヌイ試行を n 回行ったときの成功回数 X は， X ～ B(n, p)，のように２項分布に従う．X の平均と分散はそれぞれ， E[X ] = np，Var[X ] = np(1 - p)，である．
　ここで，成功確率が p₀ であるという帰無仮説，

H₀： p = p₀

の検定を考える．帰無仮説のもとでは，成功回数 X は，X ～ B(n, p₀)，と分布するので，X をその平均と標準偏差で標準化すると，中心極限定理から，

のように標準正規分布に漸近的に従う．
　これより，近似的な 5％両側検定は，標準正規分布の 97.5％分位点の z₀ = 1.96 より検定統計量 T = |z| の値が大きくなったとき帰無仮説を棄却することで得られる．なお，二項分布は離散的なので，Yates の連続性の補正（continuity correction）を行った検定統計量を用い，

のとき帰無仮説を棄却する方が近似の精度がよいと言われている．

　このように，中心極限定理を利用して，標準正規近似を行って検定を行うやり方を大標本（large sample）理論といい，コンピュータが発達する以前はもっぱら大標本理論に基づいた検定を行っていた．

比率の正規近似に基づく信頼区間

　成功確率 p のベルヌイ試行を n 回行ったとき x 回成功したとすると，成功確率は， p^{^} = x/n，と推定される．この推定値は最尤推定値である．
　成功回数 x は二項分布し，その平均は E[x ] = np，分散は Var[x ] = np(1 - p)，であるので，成功確率推定量 p^{^} の平均は E[p^{^} ] = E[x/n] = p，分散は Var[p^{^} ] = Var[x/n] = Var[x ]/n² = p(1 - p)/n，となる．これより，

と漸近的に分布するので，標準正規分布の 97.5％点の z₀ = 1.96 を用いると，近似的に

という不等式が成り立つ．これを整理すると，

という p の２次不等式を解くことに帰着する．いま，p の２次方程式の根を

とすると，この根を用い，p の 95％信頼区間は近似的に

となる．
　また，連続性の補正を行うには，成功確率の推定値 p^ を，信頼区間の下限と上限でそれぞれ

というように変えて，信頼区間が少し広くなるようにする． R では，これらの式を用いて信頼区間を構成しているようである．

　ここで，さらに近似を加えて，z₀² の項を消去すると， p の２次方程式の根は，

となるので，p の近似的な 95％信頼区間は，

と簡略化される．
　なお，この信頼区間は，成功確率推定量 p^{^} の分散において，真の成功確率 p の代わりにその推定量 p^{^} に置き換えて，Var[p^{^} ] = p^{^}(1 - p^{^})/n，とみなした場合と同じで，この信頼区間は教科書等でよく出てくる．
　簡略化された信頼区間で連続性の補正を入れるには，

として，信頼区間の幅を拡げる．

　ところで，正規近似による信頼区間の構成では，場合により信頼区間が負になったり 1 を超えることがあるが，このときは，0 と 1 で切り詰め（truncate）る．

二項確率の計算による正確な検定

　現在では二項確率が　R などのコンピュータソフトにより直接計算できるので，正規分布による近似検定を行う意味はあまりないといえる．
　いま．独立なベルヌイ試行を n 回行ったところ，成功回数が x であったとする．このとき，成功確率 p が p₀ であるという帰無仮説，とその対立仮説

H₀： p ＝ p₀， H₁： p ≠ p₀

の検定を考える．帰無仮説のもとでは成功回数の期待値は np₀ である．いま，x ＞ np₀ であるとする．このときは，まず成功回数が x 回以上である確率 Pr[X ≧ x] を計算する．これが上側確率で，

である．
　両側検定では下側確率を考える必要がある．これは，帰無仮説のもとでデータ x が得られる確率，

より，小さな確率を持つ密度を加え合わせて，

として求められる．これらから，帰無仮説の有意確率（p 値）は，

P_upper + P_lower

となり，これがが 0.05 より小さければ有意水準 5％で帰無仮説を棄却し，0.01 より小さければ有意水準 1％で帰無仮説を棄却する．
　x ＜ np₀ のときは，上の上側確率と下側確率の関係が逆になる．

二項確率のベータ分布表現による信頼区間の構成

　成功確率 p のベルヌイ試行を n 回行ったときの成功回数の確率変数を X，パラメータ，x と n - x + 1 のベータ分布に従う確率変数を Y， Y ～ beta(x，n - x + 1)，とすると，成功回数 X がある回数 x 以上になる確率は，

となる．また，成功回数が x 以下になる確率は，

となる．これらの関係は，ベータ分布の累積分布を部分積分すると，

となることを利用して求められる．

　これより，成功確率 p の 95％信頼区間は，

下限：Beta(x, n - x + 1) 分布の累積確率が 0.025 になる値，qbeta(0.025, x, n - x + 1)
上限：Beta(x + 1, n - x) 分布の累積確率が 0.975 になる値，qbeta(0.975, x + 1, n - x)

として求められる．

　ところで，F 分布とベータ分布の関係

を用いて，成功確率 p の信頼区間を F 分布を用いて構成することもできる．これは，以前は F 分布の数表が充実していたためであろう．統計ソフトが発達した現在では，ベータ分布の累積分布関数を用いて信頼区間を構成するのが普通であり，R でもベータ分布を利用している．

成功確率の正確な検定例
　成功確率の正確な検定の例として，n = 10 回の試行で x = 7 回の成功が観察されたベルヌイ試行で，

帰無仮説：H₀： p ＝ 0.4

対立仮説：H₁： p ≠ 0.4

の検定を行ってみる．
　まず上側確率を求める．帰無仮説のもとで，x = 7 回の成功が得られる確率を q₇ とすると，これ以上の成功回数 8 回，9 回，10 回が得られる確率，q₈，q₉，q₁₀，を加え合わせる．これらは，

となり，上側確率は，

となる．
　下側確率は，帰無仮説のもとで，x = 0，1，…，回の成功が得られる確率で q₇ より小さいものを， q₀，q₁，…，と求める．これらは，

となり，q₂ ＞ q₇ なので，求める下側確率は，

となる．
　これより p 値は，

P_upper + P_lower ＝ 0.05476188 + 0.0463574 ＝ 0.1011193

となる．

# 二項確率の正確な検定の R スクリプト

binom.test(7, 10, p=0.4) # 二項確率に基づく正確な検定　 up <- 1 - pbinom(6, 10, prob=0.4); up # 上側確率 pbinom(6, 10, prob=0.4, lower.tail=F) # 上側確率 lp <- pbinom(1, 10, prob=0.4); lp # 下側確率 lp + up # 両側検定 p 値 binom.test(7, 10, p=0.4, alternative="greater") # 片側検定（H1 : p > 0.4） up # 片側検定 p 値 binom.test(7, 10, p=0.4, alternative="less") # 片側検定（H1 : p < 0.4） pbinom(7, 10, prob=0.4) # 片側検定 p 値

**# 二項確率の正確な検定の R スクリプト**
binom.test(7, 10, p=0.4) # 二項確率に基づく正確な検定　 up <- 1 - pbinom(6, 10, prob=0.4); up # 上側確率 pbinom(6, 10, prob=0.4, lower.tail=F) # 上側確率 lp <- pbinom(1, 10, prob=0.4); lp # 下側確率 lp + up # 両側検定 p 値 binom.test(7, 10, p=0.4, alternative="greater") # 片側検定（H1 : p > 0.4） up # 片側検定 p 値 binom.test(7, 10, p=0.4, alternative="less") # 片側検定（H1 : p < 0.4） pbinom(7, 10, prob=0.4) # 片側検定 p 値

二項確率の近似検定と正確な検定の比較
　いま，A と B で将棋に強さに違いがあるかを調べるため，10回対戦させたところ，A の 7 勝 3 敗であった．この結果から A と B で将棋の強さに違いがあると言えるかを考えてみる．A の B に対する勝率を p とすると，帰無仮説は，A と B で将棋に強さに違いがないなら勝敗は５分５分であろうと考える．また，対立仮説は， A と B で将棋がどちらが強いか不明なので，両側とする．すなわち，

H₀： p ＝ 0.5

H₁： p ≠ 0.5

の検定を行う．ところで，比率の検定では，通常，Yates の補正を行った正規近似で検定されてきた．ここで，正確な検定との比較を行う．

**# 二項確率の検定の R スクリプト**
r <- 7; n <- 10	# 成功回数と試行回数
prop.test(r, n, p=0.5, correct=F)	# H₀：p = 0.5 の正規近似検定（補正なし）
prop.test(r, n, p=0.5)	# H₀：p = 0.5 の正規近似検定（連続性の補正）
binom.test(r, n, p=0.5)	# 二項確率に基づく正確な検定

課題：試行回数を増やして（n = 30）同様な検定を行え．

テレビ視聴率
　テレビ視聴率は，視聴率の高い番組ほど多くの視聴者が見ているので，広告宣伝の効果が高く影響力が強いと考えられている．このため，視聴率の高さが広告宣伝費用に反映されるので，テレビ会社は高い視聴率を得ようとして番組を製作している．
　ある調査会社のデータによると，関東地区では 600 世帯を対象にしているようである．
課題：ある調査会社によると，NHK 大河ドラマの関東地区世帯視聴率は24.5％であった．真の世帯視聴率の 95 ％信頼区間を求めよ．

８．適合度検定

Pearson χ² 検定

　前節の比率の検定は，χ² 分布を用いる適合度検定と大きな関係がある．ここでは n 回のベルヌイ試行で X 回成功したときに，成功確率が p₀ であるという，

帰無仮説，H₀： p ＝ p₀，対立仮説，H₁： p ≠ p₀，

の検定を考えた．そこでは，X を標準化して標準正規分布にもって行ったが，これを２乗して χ² 分布を用いることもできる．すなわち，

という関係がある．
　ところで，n 回のベルヌイ試行の結果と帰無仮説のもとでの期待値を表にすると，

	成　功	失　敗
観測度数	X	n - X
期待度数	np₀	n(1 - p₀)

となる．ここで，ピアソン（Pearson）のχ² 値，

を計算すると，

となる．つまり，χ² 値は，試行回数 n が大きくなるにつれて帰無仮説のもとで自由度 1 の χ² 分布に漸近的に従う．よって，これより検定が行える．

成功，失敗確率の適合度検定
　前節で行った A と B で将棋の強さに違いがあると言えるかの検定を，χ² 値を用いたピアソンの適合度検定で行ってみる．適合度検定では，R ではなぜか連続性の補正ができない．

# 適合度検定の R スクリプト
r <- 7; n <- 10 # 成功回数と試行回数　
x <- c(r, n-r) # 成功回数と失敗回数のベクトル　
p0 <- c(0.5, 0.5) # 成功確率と失敗確率の帰無仮説　
chisq.test(x, p=p0) # ピアソン χ² 適合度検定　

確率分布との適合度

確率分布が既知のとき

　データが想定している確率分布に適合しているかは，ピアソン（Peason）の χ² 適合度検定で行うことができる．いま，離散分布の，たとえば m = 5 のセルに対して，観測されたカウントデータと対応する想定確率が，

	セル１	セル２	セル３	セル４	セル５	計
観測度数	n₁	n₂	n₃	n₄	n₅	n
想定確率分布	p₁	p₂	p₃	p₄	p₅	1

のようになっていたとする．このとき，ピアソン（Peason）の χ² 値は，

のように近似的に自由度 m - 1 の χ² 分布に従う．これにより，データが想定確率分布に適合しているかの検定が行える．検定の帰無仮説は，

H₀：データは想定確率分布に従う．

である．　この近似は n が大きく，各セルの度数 n_i が 5 以上であることが望ましい，とされている．　一方，正規分布などの連続分布では，適当に階級分けして離散化すればこの検定が行える．ただし，階級分けは任意なので，階級分けのやり方によっては結果が異なる恐れがある．

過去の分布との比較
　以下の表のある学校での国語テストの得点分布が，過去の得点分布と同じであるかを調べたい．セルの数が 6 なので，ピアソンの χ² 適合度検定で自由度を 5 にすればよい．

得　点	30～39	40～49	50～59	60～69	70～79	80～89
人　数	1	50	123	101	24	1
過去の確率分布	0.02	0.25	0.32	0.29	0.11	0.01

# 国語得点分布適合度検定の R スクリプト
x <- c(35, 45, 55, 65, 75, 85) # 階級の代表値　
y <- c(1, 50, 123, 101, 24, 1) # 観測度数　
p0 <- c(0.02, 0.25, 0.32, 0.29, 0.11, 0.01) # 過去の確率分布　
chisq.test(y, p=p0) # χ² 適合度検定　
s <- sum(y) # 　
plot(x, y, type="h", xlab="得点", ylab="度数") # 観測度数のグラフ　
points(x, s*p0, type="b", col="red") # 期待度数のグラフ　

legend(68,113, c("観測度数", "期待度数"), lty=1, col=c("black","red"))
title(main="国語得点データと過去の得点分布との比較") # 　

**# 国語得点分布適合度検定の R スクリプト**
x <- c(35, 45, 55, 65, 75, 85)	# 階級の代表値
y <- c(1, 50, 123, 101, 24, 1)	# 観測度数
p0 <- c(0.02, 0.25, 0.32, 0.29, 0.11, 0.01)	# 過去の確率分布
chisq.test(y, p=p0)	# χ² 適合度検定
s <- sum(y)	#
plot(x, y, type="h", xlab="得点", ylab="度数")	# 観測度数のグラフ
points(x, s*p0, type="b", col="red")	# 期待度数のグラフ
legend(68,113, c("観測度数", "期待度数"), lty=1, col=c("black","red"))
title(main="国語得点データと過去の得点分布との比較")	#

二項分布
　二項分布の場合にピアソン（Peason）の χ² 値がどの程度 χ² 分布で近似されるかをみてみる．成功確率 p = 1/3 の試行を m = 6 回行ったときの成功回数 x の分布は二項分布 B(6, 1/3) に従う．これを n = 100 回行ったとすると，期待度数は，

成功回数	0 回	1 回	2 回	3 回	4 回	5 回	6 回	計
期待観測度数	8.8	26.3	32.9	22.0	8.2	1.7	0.1	100

となり，成功回数 5 回と 6 回は期待度数が小さいのでこれらをまとめ，成功回数 4 回以上とする．すると，成功回数のセル数は 5 となる．　n = 100 回を 1 回行うと観測度数が得られるので，χ² 値が計算される．これを， N = 10000 回行うと，帰無仮説が正しいときの χ² 値の分布がシミュレートされる． χ² 値のヒストグラムを自由度 5 - 1 = 4 の χ² 分布と比べればよい．　また，自由度 4 の χ² 分布の 95％分位点の値と比較して名目上の第１種の過誤が達成されているかをみる．

# 二項分布適合度検定シミュレーションの R スクリプト
N <- 10000 # シミュレーション回数　
m <- 6 # 試行１回でのベルヌイ試行の数　
n <- 100 # 二項確率の試行数　
p0 <- dbinom(x=0:6, size=6, p=1/3) # 二項確率 B(6, 1/3) 　
w <- rbinom(n, m, p=1/3) # n = 100 回の試行での成功回数（観測値）の系列）
tw <- table(factor(w,levels=0:m)) # 成功回数の表（データ）　
r <- chisq.test(tw, p=p0) # データと二項確率との適合度検定　
p1 <- c(p0[1:4], p0[5]+p0[6]+p0[7]) # 少ない確率をまとめた二項確率　
y <- NULL # χ² 値のベクトル
for(i in 1:N){ # 　
w <- rbinom(n, m, p=1/3) # n = 100 回の試行での成功回数（観測値）の系列）
fw <- table(factor(w,levels=0:m)) # 成功回数の表（データ）　
tw <- c(fw[1:4], fw[5]+fw[6]+fw[7]) # 成功回数 5，6，7 をまとめる
r <- chisq.test(tw, p=p1) # 適合度検定の χ² 値
y <- c(y, r$statistic) # 　
} # 　
my <- ceiling(max(y)) # y の最大値を超える最小の整数
mean(y) # 平均（真値 = 4）　
var(y) # 分散（真値 = 8）　

hist(y, breaks=seq(0,my,by=0.5),freq=FALSE,xlim=c(0,20),xlab="χ2 値",ylab="頻度",main="")
x <- seq(0,20, by=0.1) # 　
curve(dchisq(x, 4), 0, 20, add=T, col=2) # 自由度 4 χ² 分布　
title(main="二項分布の適合度と自由度 4 カイ自乗分布") # 　
q0 <- qchisq(0.95,df=4) # χ²(4) の 95％点　
segments(q0,0, q0,0.1, col="red") # 棄却域　
pv <- length(y[y>q0])/N # 　
s <- paste("有意水準 = ", pv) # 帰無仮説の棄却率（真値 = 0.05）　
text(15,0.06, s) # 　

確率分布のパラメータをデータから推定する場合

　確率分布のタイプ（二項分布やポアソン分布など）は想定できるが，パラメーターはデータから推定することが普通であろう．このときは，推定されたパラメーターのもとでの推定確率分布を用いて，セル数が m = 5 のときは，

	セル１	セル２	セル３	セル４	セル５	計
観測度数	n₁	n₂	n₃	n₄	n₅	n
推定確率分布	p^{^}₁	p^{^}₂	p^{^}₃	p^{^}₄	p^{^}₅	1

のような表ができる．推定したパラメーターの数が k であったとすると，このとき，ピアソン（Peason）の χ² 値は，

のように近似的に自由度 m - k - 1 の χ² 分布に従う．これにより，データが推定確率分布に適合しているかの検定が行える．

二項分布
　成功確率を推定する二項分布の場合にピアソン（Peason）の χ² 値がどの程度 χ² 分布で近似されるかをみてみる．成功確率 p = 1/3 の試行を m = 6 回行ったときの成功回数 x の分布は二項分布 B(6, 1/3) に従う．これを n = 100 回行ったとすると，たとえば，

成功回数	0 回	1 回	2 回	3 回	4 回	5 回	6 回	計
観測度数	12	18	25	32	13	0	0	100

のようなデータがえられる．このときの成功確率の推定値は， p^{^} = Σ x_i / (mn) = 0.36，となる．この推定確率から二項分布の確率分布を出し，χ² 値を計算する．この場合，成功確率をデータに合わせているので，先ほどの p = 1/3 よりは χ² 値は小さくなるはずである．このため，χ² 分布の自由度を小さくする必要がある．
　成功回数 5 回と 6 回は観測度数が小さいのでこれらをまとめ，成功回数 4 回以上とする．すると，成功回数のセル数は 5 となる．また，推定パラメーター数は 1 である．
　n = 100 回を 1 回行うと観測度数が得られるので，χ² 値が計算される．これを， N = 10000 回行うと，帰無仮説が正しいときの χ² 値の分布がシミュレートされる． χ² 値のヒストグラムを自由度 5 - 2 = 3 の χ² 分布と比べればよい．　また，自由度 3 の χ² 分布の 95％分位点の値と比較して名目上の第１種の過誤が達成されているかをみる．

# 二項分布適合度検定シミュレーションの R スクリプト
N <- 10000 # シミュレーション回数　
m <- 6 # 試行１回でのベルヌイ試行の数　
n <- 100 # 二項確率の試行数　
y <- NULL # χ² 値のベクトル
for(i in 1:N){ # 　
w <- rbinom(n, m, p=1/3) # n = 100 回の試行での成功回数（観測値）の系列）
pd <- mean(w)/m # データのもとでの推定成功確率　
p0 <- dbinom(x=0:6, size=6, p=pd) # 推定成功確率のもとでの二項確率分布
p1 <- c(p0[1:4], p0[5]+p0[6]+p0[7]) # 少ない確率をまとめた二項確率
fw <- table(factor(w,levels=0:m)) # 成功回数の表（データ）　
tw <- c(fw[1:4], fw[5]+fw[6]+fw[7]) # 成功回数 5，6，7 をまとめる
r <- chisq.test(tw, p=p1) # 適合度検定の χ² 値
y <- c(y, r$statistic) # 　
} # 　
my <- ceiling(max(y)) # y の最大値を超える最小の整数
mean(y) # 平均（真値 = 4）　
var(y) # 分散（真値 = 8）　

hist(y, breaks=seq(0,my,by=0.5),freq=FALSE,xlim=c(0,20),xlab="χ2 値",ylab="頻度",main="")
x <- seq(0,20, by=0.1) # 　
curve(dchisq(x, 3), 0, 20, add=T, col=2) # 自由度 3 χ² 分布　
title(main="二項分布の適合度と自由度 3 カイ自乗分布") # 　
q0 <- qchisq(0.95,df=3) # χ²(3) の 95％点　
segments(q0,0, q0,0.1, col="red") # 棄却域　
pv <- length(y[y>q0])/N # 　
s <- paste("有意水準 = ", pv) # 帰無仮説の棄却率（真値 = 0.05）　
text(15,0.06, s) # 　

Weldon のサイコロ実験の適合度検定
前期であつかった Weldon のサイコロ実験の適合度検定を行ってみる． 12個のサイコロを同時に投げ，5か6の目が出た個数をカウントした．

   5,6の個数       0    1    2    3    4    5

出た回数    185       1149       3265       5475       6114       5194

   5,6の個数       6    7    8    9 10 以上    合計

出た回数    3067       1331       403       105       18    26306

5,6の個数	0	1	2	3	4	5
出た回数	185	1149	3265	5475	6114	5194
5,6の個数	6	7	8	9	10 以上	合計
出た回数	3067	1331	403	105	18	26306

# サイコロ実験適合度検定の R スクリプト
x <- 0:10 # 5,6の個数　
y <- c(185,1149,3265,5475,6114,5194,3067,1331,403,105,18) # 出た回数　
m <- sum(x*y)/sum(y) # 出た回数の平均　
pd <- m/12 # 5，6 の出る確率の推定値　
p <- dbinom(x=0:12, size=12, p=1/3) # 正しいサイコロ（p = 1/3）のときの確率分布　
p0 <- c(p[1:10], p[11]+p[12]+p[13]) # 出る回数 10 回以上をまとめた確率分布　
rbind(y, p0*sum(y)) # 観測度数とモデルのもとでの期待度数　
chisq.test(y, p=p0) # データが p = 1/3 の二項分布に従っていることの検定
q <- dbinom(x=0:12, size=12, p=pd) # データから推定された確率に基づく二項分布　
q0 <- c(q[1:10], q[11]+q[12]+q[13]) # 出る回数 10 回以上にまとめた確率分布　
rbind(y, q0*sum(y)) # 観測度数とモデルのもとでの期待度数　
r <- chisq.test(y, p=q0) # データが p = pd の二項分布に従っていることの検定　
pchisq(r$statistic, df=r$parameter-1, lower.tail=F) # 有意確率（自由度１落とす）　

**# サイコロ実験適合度検定の R スクリプト**
x <- 0:10	# 5,6の個数
y <- c(185,1149,3265,5475,6114,5194,3067,1331,403,105,18)	# 出た回数
m <- sum(x*y)/sum(y)	# 出た回数の平均
pd <- m/12	# 5，6 の出る確率の推定値
p <- dbinom(x=0:12, size=12, p=1/3)	# 正しいサイコロ（p = 1/3）のときの確率分布
p0 <- c(p[1:10], p[11]+p[12]+p[13])	# 出る回数 10 回以上をまとめた確率分布
rbind(y, p0*sum(y))	# 観測度数とモデルのもとでの期待度数
chisq.test(y, p=p0)	# データが p = 1/3 の二項分布に従っていることの検定
q <- dbinom(x=0:12, size=12, p=pd)	# データから推定された確率に基づく二項分布
q0 <- c(q[1:10], q[11]+q[12]+q[13])	# 出る回数 10 回以上にまとめた確率分布
rbind(y, q0*sum(y))	# 観測度数とモデルのもとでの期待度数
r <- chisq.test(y, p=q0)	# データが p = pd の二項分布に従っていることの検定
pchisq(r$statistic, df=r$parameter-1, lower.tail=F)	# 有意確率（自由度１落とす）

12人の兄弟中の女児数のデータ
19世紀末のドイツの病院のデータによると，同じ両親で12人きょうだいがいる6155家族の女児数の数は以下のようであった．

   女児数       0    1    2    3    4    5    6

度数    7       45       181       478       829       1112       1343

   女児数       7    8    9    10    11    12    合計

度数    1033       670       286       104       24       3    6155

女児数	0	1	2	3	4	5	6
度数	7	45	181	478	829	1112	1343
女児数	7	8	9	10	11	12	合計
度数	1033	670	286	104	24	3	6155

課題：「Weldon のサイコロ実験」と同様に「12人の兄弟中の女児数のデータ」に対して二項分布に対する適合度検定を行え．

他の分布

　ポアソン分布や負の二項分布などの離散分布では，前節の二項分布のときと同様に適合度検定が行える．正規分布などの連続型分布では，多少恣意的にはなるが離散化を行えば適合度検定が行える．

虫歯数データの負の二項分布との適合度検定
　虫歯数データに対し，負の二項分布をモーメント法でパラメータ推定をしたところ， n^{^}₀ = 11.3315，p^{^}₀ = 0.7729，が得られ，最尤法によるパラメーター推定値として，n^{^} = 8.005， p^{^} = 0.7037，が得られた．これらの推定値をパラメーターとする負の二項分布の確率密度にデータ数をかけたものを期待度数とし，自由度をパラメーター分の 2 つ減らした χ² 適合度検定を行い，どちらの推定値がよく適合しているかを比較する．

小学生１人あたりの虫歯数
虫歯の数 0 1 2 3 4 5 6 7 8以上

児童の数 4 9 16 13 9 7 5 4 3

小学生１人あたりの虫歯数
虫歯の数	0	1	2	3	4	5	6	7	8以上
児童の数	4	9	16	13	9	7	5	4	3

# 虫歯の数モーメント法推定の適合度検定の R スクリプト
y <- c(4,9,16,13,9,7,5,4,3) # 児童数　
n <- 11.3315 # モーメント法による n の推定値　
p <- 0.7729 # モーメント法による p の推定値　
r <- length(y)-2 # 　
q <- dnbinom(0:r, size=n, prob=p) # 虫歯数 0 から 7 までの確率　
q[r+2] <- pnbinom(r, size=n, prob=p, lower.tail=F) # 虫歯数 8 以上の確率　
rbind(y, q*sum(y)) # 観測値と期待度数　
r <- chisq.test(y, p=q) # χ² 適合度検定　
pchisq(r$statistic, df=r$parameter-2, lower.tail=F) # 有意確率（自由度 2 落とす）　

**# 虫歯の数モーメント法推定の適合度検定の R スクリプト**
y <- c(4,9,16,13,9,7,5,4,3)	# 児童数
n <- 11.3315	# モーメント法による n の推定値
p <- 0.7729	# モーメント法による p の推定値
r <- length(y)-2	#
q <- dnbinom(0:r, size=n, prob=p)	# 虫歯数 0 から 7 までの確率
q[r+2] <- pnbinom(r, size=n, prob=p, lower.tail=F)	# 虫歯数 8 以上の確率
rbind(y, q*sum(y))	# 観測値と期待度数
r <- chisq.test(y, p=q)	# χ² 適合度検定
pchisq(r$statistic, df=r$parameter-2, lower.tail=F)	# 有意確率（自由度 2 落とす）

課題： 虫歯の数データに，負の二項分布のパラメーターの最尤推定値をあてはめたときの適合度検定を行い，モーメント法との比較を行え．

英語得点データの正規分布との適合度検定
　英語得点データの標本平均と標本分散を，正規分布の平均と分散とみなして正規分布をあてはめた． Q - Q プロットでみると，下の方の得点分布とあてはまりが悪いようにみえる．
　χ² 適合度検定を行うため，40 点以下と 41 点から 50 点，70 点以上をそれぞれ１つの階級とし， 51 点から 70 点までは 5 点きざみの階級にして，全部で 7 セルに離散化した．すると，以下の表のようにまとめられた．これに合わせて，正規分布を 7 区間に区切り，各区間ごとの確率を累積分布関数を用いて計算すれば，期待度数を求めることができる．χ² 適合度検定では，自由度をパラメーター分（平均，分散）の 2 つ減らす．

英語得点データ
   得点    40 未満 41 ～ 50 51 ～ 55 56 ～ 60 61 ～ 65 66 ～ 70 70 以上

人数    6        9        12        13        17        14        9

期待度数    4.33        13.95        11.91        13.60        12.93        10.24        13.06

英語得点データ
得点	40 未満	41 ～ 50	51 ～ 55	56 ～ 60	61 ～ 65	66 ～ 70	70 以上
人数	6	9	12	13	17	14	9
期待度数	4.33	13.95	11.91	13.60	12.93	10.24	13.06

　しかしながら，英語得点データのように正規分布などの連続分布をデータにあてはめた場合のあてはまりの良さを調べる場合，χ² 適合度検定は，離散化に恣意性が入るのであまり薦められない．連続型データの場合は，次に述べる Kolmogorov-Smirnov 検定を用いるのが普通である．

英語得点データの正規分布との Kolmogorov-Smirnov 検定
　累積分布関数 F(・) を持つと想定される母集団から，大きさ n の無作為標本，X₁，…，X_n，が抽出されたとする．このとき，

を経験（標本）累積分布関数という．この経験累積分布関数と母集団の累積分布関数との最大偏差，

を検定統計量とする．　標本が累積分布関数 F(・) を持つ母集団から抽出されたならば，最大偏差 D はそれ程大きくないと考えられる．帰無仮説，

H₀：標本は分布 F(・) から由来した

を最大偏差 D を用いて検定する方法を，Kolmogorov-Smirnov １標本検定と言う．Kolmogorov-Smirnov 検定には，２つの経験累積分布関数が同じ分布から由来したかを検定する Kolmogorov-Smirnov ２標本検定もある．

# 英語得点の正規分布適合度検定の R スクリプト
eigo <- c( 36,70,56,68,76,60,50,63,62,42,64,60,50,68,71,67, # 英語得点データ

50,65,67,57,72,64,61,66,46,80,46,51,59,32,55,65, 65,52,57,64,23,57,53,54,38,71,57,69,77,61,51,64,

63,43,65,61,51,69,72,68,53,66,68,58,73,65,62,67, 47,81,47,52,59,33,56,66,67,52,58,65,24,58,54,55)

d <- 5 # ヒストグラムの階級幅　

op <- par(mfrow = c(1, 2)) # グラフを横に２つ並べて表示　

hist(eigo, breaks=seq(0, 100, by=d), xlab="英語得点", ylab="頻度", main="")

n <- length(eigo) # データ数　
m <- mean(eigo) # 平均　
s <- sd(eigo) # 標準偏差　
x <- 0:100 # 　
curve(n*d*dnorm(x, m, s), 0, 100, add=TRUE, col="red") # 推定正規分布重ねて表示　
title(main="英語得点のヒストグラム") # 　

qqnorm(eigo, xlab="正規分布分位点", ylab="英語得点分位点", main="")

qqline(eigo, col=2) # Q-Q プロット　
title(main="正規 Q - Q プロット") # 　
par(op) # 　
a <- c(0,40,50,55,60,65,70,100) # 階級の区切り点の定義　
b <- hist(eigo, breaks=a) # 各セルに入る人数の計算　
y <- b$counts # セルの観測度数　
l <- length(a) # 　
p <- pnorm(a[2:(l-1)], mean=m, sd=s) # 階級の区切りまでの累積確率　
ps <- c(0, p) # 　
pe <- c(p, 1) # 　
q <- pe-ps # セルの確率分布　
rbind(y, sum(y)*q) # セルの観測度数と期待度数　
r <- chisq.test(y, p=q) # χ² 適合度検定　
pchisq(r$statistic, df=r$parameter-2, lower.tail=F) # 有意確率（自由度 2 落とす）　

plot(ecdf(eigo), xlab="英語得点", do.points=F, ylab="累積確率", main="")

curve(pnorm(x, m, s), 0, 100, add=TRUE, col="red") # 正規分布の累積分布　
title(main="英語得点の標本累積分布") # 　
ks.test(eigo, "pnorm", m, s) # Kolmogorov-Smirnov １標本検定　

**# 英語得点の正規分布適合度検定の R スクリプト**
eigo <- c( 36,70,56,68,76,60,50,63,62,42,64,60,50,68,71,67,	# 英語得点データ
50,65,67,57,72,64,61,66,46,80,46,51,59,32,55,65, 65,52,57,64,23,57,53,54,38,71,57,69,77,61,51,64,
63,43,65,61,51,69,72,68,53,66,68,58,73,65,62,67, 47,81,47,52,59,33,56,66,67,52,58,65,24,58,54,55)
d <- 5	# ヒストグラムの階級幅
op <- par(mfrow = c(1, 2))	# グラフを横に２つ並べて表示
hist(eigo, breaks=seq(0, 100, by=d), xlab="英語得点", ylab="頻度", main="")
n <- length(eigo)	# データ数
m <- mean(eigo)	# 平均
s <- sd(eigo)	# 標準偏差
x <- 0:100	#
curve(nddnorm(x, m, s), 0, 100, add=TRUE, col="red")	# 推定正規分布重ねて表示
title(main="英語得点のヒストグラム")	#
qqnorm(eigo, xlab="正規分布分位点", ylab="英語得点分位点", main="")
qqline(eigo, col=2)	# Q-Q プロット
title(main="正規 Q - Q プロット")	#
par(op)	#
a <- c(0,40,50,55,60,65,70,100)	# 階級の区切り点の定義
b <- hist(eigo, breaks=a)	# 各セルに入る人数の計算
y <- b$counts	# セルの観測度数
l <- length(a)	#
p <- pnorm(a[2:(l-1)], mean=m, sd=s)	# 階級の区切りまでの累積確率
ps <- c(0, p)	#
pe <- c(p, 1)	#
q <- pe-ps	# セルの確率分布
rbind(y, sum(y)*q)	# セルの観測度数と期待度数
r <- chisq.test(y, p=q)	# χ² 適合度検定
pchisq(r$statistic, df=r$parameter-2, lower.tail=F)	# 有意確率（自由度 2 落とす）
plot(ecdf(eigo), xlab="英語得点", do.points=F, ylab="累積確率", main="")
curve(pnorm(x, m, s), 0, 100, add=TRUE, col="red")	# 正規分布の累積分布
title(main="英語得点の標本累積分布")	#
ks.test(eigo, "pnorm", m, s)	# Kolmogorov-Smirnov １標本検定

課題： 英語得点データでは 25 点以下の 2 名が他の集団から離れているようにみえる．この 2 名を異常値（out-lier）としてデータから除き，正規分布との適合度検定を行い，あてはまりがどう変化したかを考察せよ．

参考文献（古い順）

Introduction to the Theory of Statistics, Mood, A. M., Graubill, F. A. & Boes, D. C., 1974, McGRAW-HILL
工学のためのデータサイエンス入門－フリーな統計環境Rを用いたデータ解析－，間瀬茂ら，2004，数理工学社
実践生物統計学－分子から生態まで－（第 1 章，第 2 章），東京大学生物測定学研究室編（大森宏ら）， 2004，朝倉書店
The R Tips データ解析環境 R の基本技・グラフィックス活用集，船尾暢男，2005，九天社
R で学ぶデータマインニング I －データ解析の視点から－，熊谷悦生・船尾暢男，2007，九天社
R で学ぶデータマインニング II －シミュレーションの視点から－，熊谷悦生・船尾暢男，2007，九天社

engei <- read.csv("engei.csv")	# csv データ読み込み
engei	# データの表示
x <- engei[,3]	# 園芸療法前の MSE 得点を x に格納（列数で指定）
x <- engei$療法前	# 園芸療法前の MSE 得点を x に格納（変数名で指定）
y <- engei$療法後	# 園芸療法後の MSE 得点を x に格納（変数名で指定）
t.test(y, x, paired=TRUE)	# １母集団 t 検定，paired = TRUE で対標本を指定
d <- y - x	# 療法後－療法前で療法の効果をみる
d	# 療法の効果の表示
t.test(d)	# １母集団 t 検定，先ほどと同じ検定
n <- length(d)	# 標本の大きさ（サンプルサイズ）
mean(d)	# 標本平均
sd(d)	# 標本標準偏差
dv <- n - 1	# 標本の自由度
t <- sqrt(n)*mean(d)/sd(d)	# 効果がないとの帰無仮説のもとでの t 値
t	# 検定統計量 t 値の表示
2*(1 - pt(t, df=dv))	# 両側検定の p 値
t0 <- qt(0.975, df=dv)	# 両側 5 ％検定の閾値
dw <- t0*sd(d)/sqrt(n)	# 95％信頼区間の幅
mean(d)-dw	# 95％信頼区間の下限
mean(d)+dw	# 95％信頼区間の上限
t.test(d, alternative="greater")	# 片側検定

母集団 A からの標本の偏差平方和：	S_A＝（n_A－1）s_A²
母集団 B からの標本の偏差平方和：	S_B＝（n_B－1）s_B²
母集団 A，B 全体での偏差平方和：	S ＝ S_A ＋ S_B ＝（n_A－1）s_A²＋（n_B－1）s_B²
母集団 A，B 共通の標本分散：

r <- 7; n <- 10	# 成功回数と試行回数
x <- c(r, n-r)	# 成功回数と失敗回数のベクトル
p0 <- c(0.5, 0.5)	# 成功確率と失敗確率の帰無仮説
chisq.test(x, p=p0)	# ピアソン χ² 適合度検定

x <- c(35, 45, 55, 65, 75, 85)	# 階級の代表値
y <- c(1, 50, 123, 101, 24, 1)	# 観測度数
p0 <- c(0.02, 0.25, 0.32, 0.29, 0.11, 0.01)	# 過去の確率分布
chisq.test(y, p=p0)	# χ² 適合度検定
s <- sum(y)	#
plot(x, y, type="h", xlab="得点", ylab="度数")	# 観測度数のグラフ
points(x, s*p0, type="b", col="red")	# 期待度数のグラフ
legend(68,113, c("観測度数", "期待度数"), lty=1, col=c("black","red"))
title(main="国語得点データと過去の得点分布との比較")	#

N <- 10000	# シミュレーション回数
m <- 6	# 試行１回でのベルヌイ試行の数
n <- 100	# 二項確率の試行数
p0 <- dbinom(x=0:6, size=6, p=1/3)	# 二項確率 B(6, 1/3)
w <- rbinom(n, m, p=1/3)	# n = 100 回の試行での成功回数（観測値）の系列）
tw <- table(factor(w,levels=0:m))	# 成功回数の表（データ）
r <- chisq.test(tw, p=p0)	# データと二項確率との適合度検定
p1 <- c(p0[1:4], p0[5]+p0[6]+p0[7])	# 少ない確率をまとめた二項確率
y <- NULL	# χ² 値のベクトル
for(i in 1:N){	#
w <- rbinom(n, m, p=1/3)	# n = 100 回の試行での成功回数（観測値）の系列）
fw <- table(factor(w,levels=0:m))	# 成功回数の表（データ）
tw <- c(fw[1:4], fw[5]+fw[6]+fw[7])	# 成功回数 5，6，7 をまとめる
r <- chisq.test(tw, p=p1)	# 適合度検定の χ² 値
y <- c(y, r$statistic)	#
}	#
my <- ceiling(max(y))	# y の最大値を超える最小の整数
mean(y)	# 平均（真値 = 4）
var(y)	# 分散（真値 = 8）
hist(y, breaks=seq(0,my,by=0.5),freq=FALSE,xlim=c(0,20),xlab="χ2 値",ylab="頻度",main="")
x <- seq(0,20, by=0.1)	#
curve(dchisq(x, 4), 0, 20, add=T, col=2)	# 自由度 4 χ² 分布
title(main="二項分布の適合度と自由度 4 カイ自乗分布")	#
q0 <- qchisq(0.95,df=4)	# χ²(4) の 95％点
segments(q0,0, q0,0.1, col="red")	# 棄却域
pv <- length(y[y>q0])/N	#
s <- paste("有意水準 = ", pv)	# 帰無仮説の棄却率（真値 = 0.05）
text(15,0.06, s)	#

x <- 0:10	# 5,6の個数
y <- c(185,1149,3265,5475,6114,5194,3067,1331,403,105,18)	# 出た回数
m <- sum(x*y)/sum(y)	# 出た回数の平均
pd <- m/12	# 5，6 の出る確率の推定値
p <- dbinom(x=0:12, size=12, p=1/3)	# 正しいサイコロ（p = 1/3）のときの確率分布
p0 <- c(p[1:10], p[11]+p[12]+p[13])	# 出る回数 10 回以上をまとめた確率分布
rbind(y, p0*sum(y))	# 観測度数とモデルのもとでの期待度数
chisq.test(y, p=p0)	# データが p = 1/3 の二項分布に従っていることの検定
q <- dbinom(x=0:12, size=12, p=pd)	# データから推定された確率に基づく二項分布
q0 <- c(q[1:10], q[11]+q[12]+q[13])	# 出る回数 10 回以上にまとめた確率分布
rbind(y, q0*sum(y))	# 観測度数とモデルのもとでの期待度数
r <- chisq.test(y, p=q0)	# データが p = pd の二項分布に従っていることの検定
pchisq(r$statistic, df=r$parameter-1, lower.tail=F)	# 有意確率（自由度１落とす）

y <- c(4,9,16,13,9,7,5,4,3)	# 児童数
n <- 11.3315	# モーメント法による n の推定値
p <- 0.7729	# モーメント法による p の推定値
r <- length(y)-2	#
q <- dnbinom(0:r, size=n, prob=p)	# 虫歯数 0 から 7 までの確率
q[r+2] <- pnbinom(r, size=n, prob=p, lower.tail=F)	# 虫歯数 8 以上の確率
rbind(y, q*sum(y))	# 観測値と期待度数
r <- chisq.test(y, p=q)	# χ² 適合度検定
pchisq(r$statistic, df=r$parameter-2, lower.tail=F)	# 有意確率（自由度 2 落とす）