統計的な研究対象となる集合全体を母集団 (population)という。北海道での「ゆめぴりか」の収量、東北地方のいもち病被害率、などである。母集団においては、その要素が従う確率分布 (母集団分布)が想定される。たとえば、収量は正規分布、被害率は二項分布などに従うと考えられる。
一つもししくは複数の母集団に対する断定や主張は統計的仮説 (statistical hypothesis) を構成し、 これには、
などが考えられる。
統計的仮説検定で用いられる仮説は、まず、帰無仮説 (null hypothesis) として与えられる。帰無仮説は棄却(reject)されることが期待される仮説であることが多い。 帰無仮説と反対の仮説を対立仮説 (alternative hypothesis) と言う。
例えば、
を考えてみる。なお、\(\mu_A\)と\(mu_B\)は母集団における平均である。 ここで、母集団AとBは異なる処理 (薬の投与など) をしているので、実験の目的は、母集団AとBの平均は異なる (処理効果がある) ということを示したい(対立仮説が正しいことを望む)のだが、まずは、「等しい(処理効果無し)」と仮定してみようという考え方で、数学の背理法と似た論理である。
なお、帰無仮説が正しいことを望む検定もある。例えば、
などである。この帰無仮説が棄却された場合、通常の統計的解析手法は厳密には適用できない場合があるので、取り扱いがより面倒になるためである。
検定統計量 (test statistic) は、帰無仮説が正しいと仮定した時に標本の値から算出され、検定に用いられる数値で、後述の、\(z\)値、\(t\)値、\(F\)値などがある。これらの統計量に対し、帰無仮説が正しいと仮定した時の分布である、標準正規分布、\(t\)分布、\(F\)分布などの帰無分布が分かっているので、帰無仮説を棄却するかどうかを判定する (両側) 5%分位点や \(p\)値が求められる。
帰無仮説が真と仮定した帰無分布のもとで、検定統計量の値以上(もしくは未満)の値 (検定統計量の値より極端な値) が得られる確率が\(p\)値である。たとえば、帰無分布として標準正規分布を用いる \(z\)検定の両側検定において、検定統計量の値である\(z\)値が \(z=2.1\)であったとすると、標準正規分布では、\(P[Z>2.1]=0.0179\) であり、両側検定なので \(P[Z<−2.1]\)となる確率も考慮して、\(p=2×0.01786=0.0357\) と計算される。これが、\(p\)値である。
z <- 2.1 # 標準正規分布が2.1以上の値をとるとき
p <- 1 - pnorm(z) # Figure 1の右側の黄色領域
2 * p # 黄色領域は左右にあるので、2倍する
## [1] 0.03572884
Figure 1: z = 2.1のときの\(p\)値に対応する面積
\(p\)値は(数学的に厳密な表現ではないが)帰無仮説が正しいとしたときに標本のようなデータが得られる確率、もしくは帰無仮説が正しい確率、と解釈するとイメージがつかみやすい。帰無仮説として「2つの母集団の平均が等しい」を考えたとき、\(p\)値が十分小さい(ほとんど起こりえない)ときは、平均が等しいと仮定したことが誤りであったと判断し、帰無仮説を棄却する。このとき、2つの母集団平均には有意な差があると結論づける。 この確率がそれほど小さくない場合は、帰無仮説のもとでも計算された統計量の値が得られる場合があると考え、帰無仮説を採択し、平均が等しいと考えても良いとする。ただし、平均は等しかった、と言い切ることはできず、平均には有意な差が認められなかった、という結論になる。
棄却か採択かの判断の基準となる確率を有意水準 (significance level) といい、切りの良い数値として 5%や 1%が良く用いられている。帰無分布の5%や1%などの分位点があらかじめ与えられていれば、検定統計量の値を分位点と比較することで、帰無仮説を棄却するかどうかを決定することができる。\(z\)検定においては、標準正規分布の両側\(5%\) 点が 1.960なので、有意水準5%の検定では、検定統計量である \(z\)値が以下の領域\(R\)に含まれたとき、 \[ R = \{z < -1.96 \} \cup \{z \ge 1.96\} \] すなわち、\(z∈R\)であるとき帰無仮説は棄却される。この領域 \(R\)を棄却域(rejection region) という。また、棄却域\(R\)の補集合 \[ Rc={−1.96\le z<1.96} \] を採択域 (acceptance region) と言う。
なお、先ほどの例では、\(z=2.1\)であり、標準正規分布の両側1%点は2.576、5%点は1.960なので、1% 有意ではないが、5%の棄却域に含まれるので、5%有意であると結論づけられる。なお、\(p\) 値が0.0357だったことからも同様の結論を導くことができる。
母集団母数(パラメータ)の検定では、薬投与などの処理を行った集団(処理群)平均\(\mu_A\) が,薬を投与しない集団(対照群)の平均\(\mu_B\)より小さくなることはないことが事前に わかっているような場合がある。このようなとき、 帰無仮説、\(H_0:\mu_A=\mu_B\)、対立仮説、\(H_1:\mu_A>\mu_B\)とすることができる。これを片側検定という。これは、事前情報より、\(\mu_A<\mu_B\)となる可能性をまったく考えない場合である。このため検定には、片側5%点や1%点を用いる。また、\(\mu_A>\mu_B\)となる可能性をまったく考えない 帰無仮説、\(H_0:\mu_A=\mu_B\)、対立仮説、\(H_1:\mu_A<\mu_B\)となる片側検定を行う場合もある。
一般には、事前情報が無いと考え、両側検定のもとでの 帰無仮説、\(H_0:\mu_A=\mu_B\)、対立仮説、\(H_1:\mu_A\ne\mu_B\) で検定を行う。このときは、両側5%点や1%点を用いる。両側検定の方が片側検定より帰無仮説を棄却しにくい (有意になりにくい)。
統計的仮説検定は、帰無仮説を棄却するか採択するかのいずれかであるが、 統計量は分布をもつので、この判定には間違いが起こることがある。 以下のように、この過誤には 2 種類がある。
仮説の棄却 (Reject) | 仮説の受託 (Accept) | |
---|---|---|
仮説が真 (True) のとき | 第1種の過誤 (Type I error) 確率 \(\alpha\) |
正解 |
仮説が偽 (False) のとき | 正解 検出力 \(1-\beta\) |
第2種の過誤 (Type II error) 確率 \(\beta\) |
第1種の過誤の確率\(\alpha\)が有意水準である。また、第2種の過誤の確率を\(\beta\)としたとき、 仮説が偽のとき正しく仮説を棄却する確率\(1−\beta\)を検出力、もしくは検定力(power)という。 良い検定は,第1種の過誤を固定したもとで検出力の高い検定方式である。
第1種の過誤の確率(赤)、第2種の過誤の確率(青)、検出力(緑)に対応する面積
母集団が従うと想定した確率分布を \(f(x|\theta)\) としたとき、\(\theta\) を分布パラメータ (母数) という。正規分布なら \(\theta = (\mu, \sigma^2)\) であり、二項分布なら \(\theta = (n, p)\) である。母集団の分布パラメータ \(\theta\) を母集団からの大きさ\(n\) の無作為標本(random sample)\(X_1, \ldots,X_n\)、を用いて推定することを点推定 (point estimation) と言う。
点推定を行う方法として、モーメント法(moment matching)と最尤法(maximum likelihood)が知られている。
モーメント法は、推定量を計算する最も自然なものである。これは、分布の \(k\) 次モーメントと標本の \(k\) 次モーメント \[ \bar{X^k} = \frac{1}{n} X_i^k \] を等しいとおくことで \(\theta\) の推定量(estimator)を求める方法である。なお、標本の実現値\(x_1,\ldots,x_n\) が得られたときは \(\theta\) の推定値(estimate)と言い、両者は区別される。
ロシア人の統計学者Bortkewitchが、プロイセン(ドイツ)陸軍の10の騎兵連隊の中で、1875年から1894年にかけての20年間におけるのべ200連隊において、1年間で「馬に蹴られて死んだ兵士の数」を調査したところ、以下の表のデータを得た。
死亡者数 | 0 | 1 | 2 | 3 | 4 | 5以上 | 計 |
---|---|---|---|---|---|---|---|
連隊数 | 109 | 65 | 22 | 3 | 1 | 0 | 200 |
このデータ分布の平均と分散 \[ \mu = {\rm E}[X] = \sum_i x_i p_i\\ \sigma^2 = {\rm Var}[X] = \sum_i (x_i-\mu)^2 p_i \] を計算し、これがポアソン分布 \[ p(x|\lambda) = \frac{\lambda^x}{x!} e^{-\lambda}, \ x =0, 1,2, \cdots, \ \mu=\sigma^2 = \lambda \] に当てはまることを示し、モーメント法でパラメータを推定してみよう。
x <- 0:5 # 死亡者数
y <- c(109, 65, 22, 3, 1, 0) # 連隊数データ
n <- sum(y); n
## [1] 200
p <- y / n # データ確率分布
m <- sum(x * p); m # データ平均
## [1] 0.61
v <- sum((x - m)^2 * p); v # データ分散
## [1] 0.6079
データ分布平均と分散がほぼ等しいので、このデータ分布はポアソン分布に従っていると考えられる。ポアソン分布パラメータ \(\lambda\) のモーメント法による推定値は、データ平均を用いる場合とデータ分散を用いる場合の2通りが考えられる。すなわち、\(\lambda_1 = 0.61、\lambda_2 = 0.6079\) の2通りである。
最尤推定(maximum likelihood estimation)は、1920年代に R. A. Fisher によって提案され、現在においても、統計的推定論において最もよく用いられている手法の一つである。この方法は、データ生成モデルを母数の関数とみなしたときの、極値を探索する最適化問題として定式化できる。
いま、\(X_1,\ldots,X_n\) を \(f(x|\theta)\) の分布をもつ母集団から得られた標本だとする。標本の要素は互いに独立なので、その同時分布は、個々の確率密度関数の積で表され、 \[ f(X_1,\ldots,X_n|\theta) = \prod_i f(X_i|\theta) \] となる。この同時分布を母数 \(\theta\) の関数としてみたものを尤度(likelihood)、 \[ L(\theta|X_1,\ldots,X_n ) = \prod_i f(X_i|\theta) \] と呼ぶ。尤度 \(L(\theta|X_1,\ldots,X_n)\) を最大化する統計量 \(\hat{\theta}\) を最尤推定量(maximum likelihood estimator : MLE)と言う。
多くの場合、尤度の対数である対数尤度(log likelihood)、 \[ l(\theta| X_1,\ldots,X_n) = \log L(\theta| X_1,\ldots,X_n ) = \sum_i \log f(X_i|\theta) \] を最大化する方が容易である。これは尤度における掛け算が対数尤度では足し算になるからである。対数は単調増加関数なので、尤度と対数尤度の最適化は同じ最尤推定量になる。
具体的には、 \[ \hat{\theta} = {\rm argmax}_\theta l(\theta| X_1,\ldots,X_n) \] であり、これは、 \[ \frac{\partial^2 l(\theta| X_1,\ldots,X_n)}{\partial \theta^2} < 0 \] (対数尤度の2次微分が負になる、すなわち上に凸)の条件のもとで、 \[ \frac{\partial l(\theta| X_1,\ldots,X_n)}{\partial \theta} = 0 \] (対数尤度の1次微分が0になる、すなわち傾き0)の解として得ることができる。
ポアソン分布からの標本 \(X_1,\ldots,X_n\) に対し、尤度は、
\[ L(\lambda | X_1,\ldots,X_n) = \prod_i e^{-\lambda}\frac{\lambda^{X_i}}{X_i !} = e^{-n\lambda}\frac{\lambda^{\sum_i X_i}}{\prod_i X_i !} \] と書ける。尤度の最大を直接求めるより、対数尤度、 \[ l(\lambda | X_1,\ldots,X_n) = -n \lambda + \sum_i X_i \log\lambda - \log \prod_i X_i ! \] を考えた方が計算が容易になる。上式を\(\lambda\) で微分して、0とおくと、 \[ \frac{d l}{d \lambda} = -n + \frac{\sum_i X_i}{\lambda} = 0 \] となり、最尤推定量は、 \[ \hat{\lambda} = \frac{1}{n}\sum_i X_i = \bar{X} \] となる。これは1次のモーメント推定量と一致している。なお、2次微分は、 \[ \frac{d^2 l}{d \lambda^2} = -\frac{\sum_i X_i}{\lambda^2} < 0 \] となるので常に負(\(\because X_i\)は自然数)であり、対数尤度は上に凸であり、最尤推定量は対数尤度を最大にすることがわかる。
正規分布 \(N(\mu, \sigma^2)\) からの標本 \(X_1,\ldots,X_n\) に対し、尤度及び対数尤度は、 \[ L(\mu,\sigma^2|X_1,\ldots,X_n) = \prod_i \frac{1}{\sqrt{2 \pi \sigma^2}} e^{-\frac{(X_i-\mu)^2}{2\sigma^2}} = \Bigl(\frac{1}{2 \pi \sigma^2} \Bigr)^{n/2} e^{-\frac{1}{2 \sigma^2}\sum_i(X_i-\mu)^2} \\ l(\mu,\sigma^2|X_1,\ldots,X_n) = -\frac{n}{2}\log 2 \pi - \frac{n}{2}\log \sigma^2 - \frac{1}{2 \sigma^2} \sum_i(X_i - \mu)^2 \] となる。対数尤度を \(\mu、\sigma^2\) で偏微分して0とおくと、 \[ \frac{\partial l}{\partial \mu} = \frac{1}{\sigma^2}\sum_i(X_i - \mu) = 0 \\ \frac{\partial l}{\partial \sigma^2} = -\frac{n}{2\sigma^2}+\frac{1}{2(\sigma^2)^2} \sum_i (X_i-\mu)^2 = 0 \] となり、これらを解くと、最尤推定量 \[ \hat{\mu} = \frac{1}{n}\sum_i X_i = \bar{X}\\ \ \hat{\sigma^2} = \frac{1}{n} \sum_i (X_i-\hat{\mu})^2 = \frac{1}{n} \sum_i (X_i - \bar{X})^2 \] が得られる。
\(X_1,\ldots,X_n\) が \(f(x|\theta)\) の分布をもつ母集団から得られた標本であり、\(\hat{\theta}_n = g (X_1,\ldots,X_n)\) が母数 \(\theta\) を推定するための統計量だとする。推定量 \(\hat{\theta}_n\) は、標本の関数なので、確率変数である。確率変数であるならば、推定量はその期待値(平均)や分散、および、標本分布(sampling distribution)と呼ばれる分布を持つ。
標本分布に対する期待値が、母集団の真の値に一致する、すなわち、 \[{\rm E}[\hat{\theta}_n] = \theta\] のとき、推定量 \(\hat{\theta}_n\) は不偏(unbiased)と言う。また、 \[ b(\hat{\theta}_n ) = {\rm E}[\hat{\theta}_n ] - \theta \] を \(\hat{\theta}_n\) の偏り、バイアス(bias)と呼ぶ。
母集団から大きさ \(n\) の標本 \(X_1, \cdots, X_n\) を抽出して分布パラメータ \(\theta\) の推定量 \(\hat{\theta_n}\) 構成したとする。このとき、標本の大きさ \(n\) を大きくしていけば母数 \(\theta\) を正しく推定できることが必要である。つまり、任意の \(\epsilon\) に対し、 \[ \lim_{n \to \infty} {\bf P} [ |\hat{\theta_n}-\theta | < \epsilon ] = 1 \] が成り立つことが必要である。このような性質を持つ推定量を一致推定量 (consistent estimator) という。
不偏推定量の分散の挙動を調べることで、その一致性が確認できる。すなわち、不偏推定量 \(\hat{\theta}_n\) において、標本の大きさ \(n \rightarrow \infty\) のとき、Var\([\hat{\theta}_n] \rightarrow 0\) となるとき、その推定量は一致性があるという。
母平均 \(\mu\) の推定量として標本平均 \(\hat{\mu} = \bar{X}\) は自然な選択である。推定量 \(\bar{X}\) は、多くの分布において母平均 \(\mu\) の「最適な」推定量である。
推定量 \(\bar{X}\) は標本ごとに変化する。より正確に言うと、 \(\bar{X}\) は、観察値 \(X_i\) の分布に依存して決まる分布を持つ確率変数である。以下は、E\([X_i] = \mu\)、Var\([X_i] = \sigma^2\) が存在する限り、母集団のいかなる分布についても成り立つ。これは、\(\bar{X}\) が不偏で一致性をもつ \(\mu\) の推定量であることを示している。 \[ {\rm E}[\bar{X}] = \mu\\ {\rm Var}[\bar{X}] = \frac{\sigma^2}{n} \] 証明は以下の通りである。標本 \(X_1,\ldots,X_n\) は互いに独立で同一の分布 (母集団分布) に従う (independent and identically distributed. iid) ので、 \[ {\rm E}[\bar{X}] = {\rm E}[\frac{1}{n}\sum_i X_i] = \frac{1}{n}\sum_i {\rm E}[X_i] = \frac{1}{n}n\mu = \mu \\ {\rm Var}[\bar{X}] = {\rm Var}[\frac{1}{n}\sum_i X_i] =\frac{1}{n^2}\sum_i{\rm Var}[X_i]= \frac{1}{n^2}n\sigma^2 = \frac{\sigma^2}{n} \] である。
なお、標本の正規性\(X_i \sim N(\mu, \sigma^2)\)を仮定すると、\(\bar{X}\) の標本分布は、正規分布の再生性から、 \[ \bar{X} \sim N \Bigl(\mu, \frac{\sigma^2}{n} \Bigr) \] の正規分布に従う。
標本 \(X_1,\ldots,X_n\) に対し、母平均 \(\mu\) が既知の場合、\(\sigma^2\) の推定量は、 \[ \hat{\sigma^2} = \frac{1}{n}\sum_i(X_i - \mu)^2 \] となる。\(\mu\) が未知で、これを \(\bar{X}\) で推定する場合、 \[ \hat{\sigma^2} = s^2 = \frac{1}{n-1}\sum_i(X_i - \bar{X})^2 \] で推定される。\(s^2\) は標本分散とも呼ばれ、\(n\) でなく \(n-1\) で割るのは \(\sigma^2\) の不偏推定量にするためである。実際、 \[ {\rm E}[s^2] = {\rm E} \Bigl[\frac{1}{n-1}\sum_i(X_i - \bar{X})^2 \Bigr] = \frac{1}{n-1}{\rm E}[\sum_i(X_i - \mu)^2 -n(\bar{X}-\mu)^2 ] \\ = \frac{1}{n-1} \left\{ \sum_i {\rm E[(X_i - \mu)^2]} -n{\rm E}[(\bar{X}-\mu)^2] \right\} = \frac{1}{n-1} \{\sum_i{\rm Var}[X_i] - n{\rm Var}[\bar{X}] \} \\ = \frac{1}{n-1} \Bigl(n\sigma^2 - n\frac{\sigma^2}{n} \Bigr) = \sigma^2 \] となるので、\(s^2\) が不偏であることが分かる。
分散 \(\sigma^2\) の最尤推定量(\(n\) で割る)は不偏ではなく、\(n-1\) で割る標本分散は不偏であることをシミュレーションで確認してみよう。いま、平均 50、分散 100の正規母集団を考える。母集団平均と分散が未知であるときに、ここから大きさ \(n = 20\) の無作為標本を抽出して母分散を推定した。\(N = 10000\) 回のシミュレーションで標本分散を計算したところ、その平均値が 99.72となり、ほぼ 100なので、標本分散が不偏推定量であることが示された。一方、\(n\) で割る最尤推定量は不偏ではなかった。
n <- 20; N <- 10000 # 標本の大きさとシミュレーション回数
v1 <- v2 <- rep(NA, N)
for(i in 1:N){
x <- rnorm(n, mean=50, sd=10) # N(50, 10^2) からn個の正規乱数抽出
v1[i] <- var(x) # 標本分散
v2[i] <- var(x) * (n-1) / n # nで割る分最尤推定量
}
mean(v1)
## [1] 99.76006
mean(v2)
## [1] 94.77206
標本分散の一致性についてもシミュレーションで確認してみる。一致性は、1回の推定に用いる標本の大きさ\(n\)を無限大にしたときに、母集団のパラメータの値(この場合は、母分散)に一致するというものである。そこで、\(n\)を大きくしながら、推定される標本分散の挙動を確認してみる。
n <- 1000
v1 <- rep(NA, n)
for(i in 2:n) {
x <- rnorm(i, mean=50, sd=10)
v1[i] <- var(x)
}
plot(v1, type = "l")
標本の大きさ\(n\)が大きくなるに連れ、母分散(100)に近づいていくことがわかる。
点推定量 \(\hat{\theta}\) の標本分布が連続的な場合、必ず \({\bf P}[\hat{\theta}=\theta] = 0\)、となる。つまり、推定量が未知母数に完全に一致する確率は 0 である。そこで、2 つの推定量 \(L = L(X_1,\ldots,X_n)\) と \(U = U(X_1,\ldots,X_n)\) を考え、\([L, U]\) の区間に \(\theta\) が \(1-\alpha\) の確率で含まれるようにすることが考えられる。このとき、区間 \([L, U ]\) は \(\theta\) の \((1-\alpha) \times 100\%\) 信頼区間とよばれる。例えば、\(\alpha= 0.05\) のときは、95% 信頼区間が構成される。
まず、正規分布に従う母集団の母平均 \(\mu\) に対する検定を考える。過去のデータなどから母分散 \(\sigma^2\) の値が既知である場合と、標本から母分散を推定する場合とで検定方式は異なる。これは、母分散を標本から推定すると、その推定誤差が検定方式に影響を与えるからである。
分散既知の正規分布 \(N(\mu, \sigma^2)\) から大きさ \(n\) 標本、\(X_1, \cdots, X_n\) を抽出したとき、
\[ {\rm H}_0 : \mu = \mu_0\\ {\rm H}_1 : \mu \ne \mu_0 \] の両側検定問題を考えてみよう。母平均 \(\mu\) の推定量は標本平均 \(\bar{X}\) が最適であった。また、その標本分布は、\(\bar{X} \sim N(\mu, \sigma^2/n)\) であり、これを標準化した変量 \(Z\) は、 \[ Z = \frac{\bar{X}-\mu}{\sigma/\sqrt{n}} \sim N(0,1) \] のように標準正規分布に従っている。いま、帰無仮説が正しいと仮定すると、\(\mu = \mu_0\) なので、 \[ z = \frac{\bar{X}-\mu_0}{\sigma/\sqrt{n}} \] である \(z\) 値を検定統計量として、標準正規分布の分位点と比較することで母平均に対する検定が行える。両側検定であれば、標準正規分布の \(1-\alpha/2\) 分位点を \(z_{1-\alpha/2}\) とすると、\(|z| > z_{1-\alpha/2}\) であれば、帰無仮説を有意水準 \(100\alpha \%\) で棄却する。例えば、\(\alpha=0.05\) ならば標準正規分布の \(97.5\%\) 分位点が1.96なので、両側検定の場合、\(|z| > 1.96\) であれば、帰無仮説を有意水準 \(5\%\) で棄却する。このような検定方式を \(z\) 検定という。
分散未知の正規分布 \(N(\mu, \sigma^2)\) から大きさ \(n\) 標本、\(X_1, \cdots, X_n\) を抽出したとき、
\[ 帰無仮説、{\rm H}_0 : \mu = \mu_0、 対立仮説、{\rm H}_1 : \mu \ne \mu_0 \]
の検定問題を考えてみよう。母平均 \(\mu\) の推定量は標本平均 \(\bar{X}\)、母分散 \(\sigma^2\) の推定量は標本分散 \(s^2 = \frac{1}{n-1}\sum_i (X_i - \bar{X})^2\) で推定される。このとき、 \[ t = \frac{\bar{X}-\mu}{s/\sqrt{n}} \sim t_{n-1} \] は自由度 \(n-1\) の \(t\) 分布に従う。よって、帰無仮説が真ならば、\(\mu =\mu_0\) であるので、 \[ t = \frac{\bar{X}-\mu_0}{s/\sqrt{n}} \] である \(t\) 値を検定統計量として、自由度 \(n-1\) の \(t\) 分布の分位点と比較することで母平均に対する検定が行える。両側検定であれば、 自由度 \(n-1\) の \(t\) 分布の \(1-\alpha/2\) 分位点を \(t_{n-1, 1-\alpha/2}\) とすると、\(|t| > t_{n-1, 1-\alpha/2}\) であれば、帰無仮説を有意水準 \(100\alpha \%\) で棄却する。このような検定方式を \(t\) 検定という。
\(X_1,\ldots,X_n \sim N(\mu, \sigma^2)\) であり、母分散 \(\sigma^2\) が既知であるとする。標本平均を標準化して、標準正規分布に従う確率変数 \(Z\) を求めると、 \[ X_i \sim N(\mu, \sigma^2) \rightarrow \bar{X} \sim N \Bigl(\mu, \frac{\sigma^2}{n} \Bigr), \ Z = \frac{\bar{X}-\mu}{\sigma/\sqrt{n}} \sim N(0,1) \] となる。いま、標準正規分布 \(N(0, 1)\) の \(1-\alpha/2\) 分位点を \(z_{1-\alpha/2}\) とすると \[ {\bf P}[-z_{1-\alpha/2} < Z < z_{1-\alpha/2}] = 1-\alpha \\ {\bf P}[-z_{1-\alpha/2} < \frac{\bar{X}-\mu}{\sigma/\sqrt{n}} < z_{1-\alpha/2}] = 1-\alpha \\ {\bf P}[-z_{1-\alpha/2}\frac{\sigma}{\sqrt{n}} + \mu < \bar{X} < z_{1-\alpha/2}\frac{\sigma}{\sqrt{n}} + \mu] = 1-\alpha \\ {\bf P}[ \bar{X} -z_{1-\alpha/2}\frac{\sigma}{\sqrt{n}} < \mu < \bar{X} +z_{1-\alpha/2}\frac{\sigma}{\sqrt{n}} ] = 1-\alpha \] となる。最後の式が、母平均 \(μ\) の \((1-\alpha)\times 100 \%\) 信頼区間である。
\(\sigma^2\) が未知の場合、信頼区間の式の \(\sigma\) は標本標準偏差 \(s\) で置き換え、\(z\) の分位点は自由度 \(n-1\) の \(t\) 分布の分位点 \(t_{n-1, 1-\alpha/2}\) を用いる必要がある。よって、 \[ \bar{X} -t_{n-1,1-\alpha/2}\frac{s}{\sqrt{n}} < \mu < \bar{X} + t_{n-1, 1-\alpha/2}\frac{s}{\sqrt{n}} \] が母平均 \(μ\) の \((1-\alpha)\times 100\%\) 信頼区間となる。
母平均に対する両側検定は,母平均に対する信頼区間と大きな関係がある。いま、帰無仮説(\({\rm H}_0\))と対立仮説(\({\rm H}_1\))が、
\[ {\rm H}_0 : \mu = \mu_0\\ {\rm H}_1 : \mu \ne \mu_0 \] であるとする。このとき、
が成り立つ。
95%信頼区間の 95%の意味は、未知の母平均 \(\mu\) が 95%の確率でその区間に含まれる、という意味ではない。正確には以下のようになる。母集団からの標本抽出を多数回行ったとする。すると、標本ごとに母平均 \(\mu\) の推定値である標本平均と \(\mu\) の 95%信頼区間が計算される。この構成された多くの信頼区間の 95%が母平均 \(\mu\) を含む、という意味である。
95%信頼区間の意味をシミュレーションで確認してみよう。正規母集団 \(N(10, 22)\) から \(n =30\) の標本を抽出して95%信頼区間を計算する。これを \(N =10000\) 回行い、どれだけの割合の信頼区間が真の母平均 \(\mu= 10\) を含むかを調べる。
N <- 10000 # 標本抽出(サンンプリング)回数
n <- 30 # 標本の大きさ(サンプルサイズ)
mu <- 10; sigma <- 2 # 母集団の母平均と母標準偏差
alpha <- 0.05 # (1-α)100%信頼区間
tquantile <- qt(1 - alpha/2, df = n-1) # t分布分位点
tquantile # 2.04523 > 1.96,正規分布より信頼区間の幅が広くなる。
## [1] 2.04523
# 正規母集団からの標本抽出のシミュレーション
X <- matrix(rnorm(N*n, mu, sd=sigma), nrow=N, byrow=T)
xbar <- apply(X, 1, mean) # 標本抽出ごとの標本平均
s <- apply(X, 1, sd) # 標本抽出ごとの標本標準偏差
LB <- xbar - tquantile*s/sqrt(n) # 標本抽出ごとの信頼区間下限
UB <- xbar + tquantile*s/sqrt(n) # 標本抽出ごとの信頼区間上限
cover <- LB < mu & UB > mu # 信頼区間に母平均μが含まれるかの真(T)偽(F)
head(cover) # 含まれていればTRUE、含まれていなければFALSE
## [1] TRUE TRUE TRUE TRUE TRUE TRUE
sum(cover)/N # 計算された信頼区間がμを含んだ割合
## [1] 0.9491
# 信頼区間100個の表示
id <- 1:100
col <- rep("black", length(id))
col[!cover[id]] <- "red" # 信頼区間がμを含まない時は赤
matplot(cbind(LB, UB)[id,], type="n")
segments(id, LB[id], id, UB[id], col=col)
points(id, xbar[id], pch=4, col="blue", cex=0.8)
abline(h=10, col="blue")
title(main="100 confidence intervals\nRed ones failed to contain the true mean")
ある栽培条件の下で栽培されたリンゴの収穫適期における果実の重さが350gであるかどうかを調べたい。なお、その栽培条件で栽培されたリンゴの果実重の標準偏差は50であることは知られているとする。いま、25個の果実の平均の果実重 \(\bar{X}\) を用いた両側検定を有意水準5%で行うとする。
\(z_{0.975}\)を標準正規分布の97.5%分位点(両側5%点)とすると、帰無仮説\(H_0:\mu = \mu_0\)が採択されるのは、\(\bar{X}\)が以下の不等式を満たすときである。
\[ - z_{0.975} < \frac{\bar{X} - \mu_0}{\sigma/\sqrt{n}} < z_{0.975} \\ \mu_0 - z_{0.975} \frac{\sigma}{\sqrt{n}} < \bar{X} <\mu_0 + z_{0.975} \frac{\sigma}{\sqrt{n}} \]
まず、\(z_{0.975}\)を求める。
z0 <- qnorm(0.975)
z0
## [1] 1.959964
次に、\(\mu_0 = 350, \sigma = 50, n = 25\)であるので、これらを代入して採択域\(z_{0.975} \frac{\sigma}{\sqrt{n}}\)の幅(片側)計算し、それを\(\mu_0\)から減ずる・加える。
mu0 <- 350
sigma <- 50
n <- 25
d <- z0 * sigma / sqrt(n)
c(mu0 -d, mu0 + d)
## [1] 330.4004 369.5996
帰無仮説 (\(\mu_0 = 350\)) の元では、\(\bar{X} \sim N (\mu_0, \sigma^2/n)\) と分布するはずだが、実際は、\(\bar{X} \sim N (\mu_1, \sigma^2/n), \mu_1=380\) と分布する。両側 5%で帰無仮説が採択されるには、\(z_0 = (\bar{X}-\mu_0)/(\sigma/\sqrt{n}) < z_{0.975}\)、すなわち、\(\bar{X} < \mu_0 + z_{0.975} \sigma/\sqrt{n}\) となる必要がある。
この両辺から \(\mu_1\) を引き、それを\(\sigma/\sqrt{n}\) で割ると、\(z_1 = (\bar{X}-\mu_1)/(\sigma/\sqrt{n}) < (\mu_0 - \mu_1)/(\sigma/\sqrt{n}) + z_{0.975}\) となる。\(z_1\) が真の値なので、これより帰無仮説が採択される確率が求められる。
mu1 <- 380
z1 <- (mu0 - mu1) / (sigma / sqrt(n)) + z0
z1
## [1] -1.040036
標準正規分布が\(z_1\)の値をとったときの累積確率を求める。
pnorm(z1)
## [1] 0.1491616
このとき、検出力(\(1-\beta\))は、
1 - pnorm(z1)
## [1] 0.8508384
なお、Rの関数を用いると\(z_0, z_1\)を計算せずに計算をすることもできる。
例えば、問題1は、
qnorm(c(0.025, 0.975), mean = 350, sd = 50 / sqrt(25))
## [1] 330.4004 369.5996
とすれば求められる。
また、問題2は、
ub <- qnorm(0.975, mean = 350, sd = 50 / sqrt(25))
pnorm(ub, mean = 380, sd = 50 / sqrt(25))
## [1] 0.1491616
とすれば求められる。
最後にグラフを描く。
x <- seq(300, 450, 1)
d0 <- dnorm(x, mean = 350, sd = 50 / sqrt(25))
d1 <- dnorm(x, mean = 380, sd = 50 / sqrt(25))
plot(x, d0, type = "l")
lines(x, d1)
abline(h = 0)
abline(v = ub)
上の図の、どこの領域が第1の過誤や第2の過誤、検出力にあたる部分か考えてみよう。
先の計算では、両側 5%で帰無仮説が採択される場合の上限のみを計算し、第2種の過誤が生じる確率を計算していた。
ub <- qnorm(0.975, mean = 350, sd = 50 / sqrt(25))
pu <- pnorm(ub, mean = 380, sd = 50 / sqrt(25))
pu
## [1] 0.1491616
しかし、厳密には、下限よりも小さな値が得られて、帰無仮説が棄却される可能性もある。この確率は、
lb <- qnorm(0.025, mean = 350, sd = 50 / sqrt(25))
pl <- pnorm(lb, mean = 380, sd = 50 / sqrt(25))
pl
## [1] 3.525313e-07
したがって、厳密には第2種の過誤は、
pu - pl
## [1] 0.1491612
となる。plの値が小さいために、ほとんど結果は変わらないが、小数点以下7桁目は変化する。
このように、真の分布が帰無分布と離れているときには、下限を下回る確率を考慮しなくても、大きな影響はないが、真の分布が帰無分布に近くなればなるほど、下限を下回る確率を考慮する必要がある。例えば、リンゴの果実重の真の平均は360gだったとすると、
ub <- qnorm(0.975, mean = 350, sd = 50 / sqrt(25))
pu <- pnorm(ub, mean = 360, sd = 50 / sqrt(25))
pu
## [1] 0.8314633
lb <- qnorm(0.025, mean = 350, sd = 50 / sqrt(25))
pl <- pnorm(lb, mean = 360, sd = 50 / sqrt(25))
pl
## [1] 0.001538375
pu - pl
## [1] 0.829925
真の平均ga352gだったとすると、
ub <- qnorm(0.975, mean = 350, sd = 50 / sqrt(25))
pu <- pnorm(ub, mean = 352, sd = 50 / sqrt(25))
pu
## [1] 0.960793
lb <- qnorm(0.025, mean = 350, sd = 50 / sqrt(25))
pl <- pnorm(lb, mean = 352, sd = 50 / sqrt(25))
pl
## [1] 0.01538773
pu - pl
## [1] 0.9454053
すでに点推定のところで述べたように、正規母集団における母分散 \(\sigma^2\) は、母平均 \(\mu\) が未知の場合、標本分散 \(s^2\)、 \[ \hat{\sigma^2} = s^2 = \frac{1}{n-1}\sum_i(X_i - \bar{X})^2 \] で推定される。ところで、標本分散の標本分布は、 \[ \frac{(n-1)s^2}{\sigma^2} \sim \chi^2_{n-1} \] のように自由度 \(n-1\) の \(\chi^2\) 分布に従うので、これを用いて母分散の検定を構成することができる。いま、母分散に対する両側検定 \[ {\rm H}_0 : \sigma^2 = \sigma^2_0\\ {\rm H}_1 : \sigma^2 \ne \sigma^2_0 \] を考えると、帰無仮説が正しいと仮定するので、検定統計量は、\(\sigma^2\) に \(\sigma^2_0\) を代入した、 \[ \chi^2 = \frac{(n-1)s^2}{\sigma^2_0} \] として、自由度 \(n-1\) の \(\chi^2\) 分布の分位点を用いて構成できる。有意水準 \(\alpha\) の両側検定は、自由度 \(n-1\) の \(\chi^2\) 分布の \(\alpha\) 分位点を \(\chi^2_{n-1, \alpha}\) とすると、\(\chi^2\) 分布は左右対称でないので、\(\chi^2_1 = \chi^2_{n-1, \alpha/2}、\chi^2_2 = \chi^2_{n-1, 1-\alpha/2}\) としたとき、\(\chi^2 < \chi^2_1\) もしくは \(\chi^2_2 < \chi^2\) が成り立つとき帰無仮説を棄却する。
母分散に対する検定で示したように、標本分散 \(s^2\) の標本分布は、自由度 \(n-1\) の \(\chi^2\) 分布 に従うので、これを用いて母分散の信頼区間を構成することができる。すなわち、 \[ {\bf P}[\chi^2_1 < \chi^2_{n-1} < \chi^2_2] = 1-\alpha \\ {\bf P}[\chi^2_1 < \frac{(n-1)s^2}{\sigma^2} < \chi^2_2] = 1-\alpha \] となる。これを整理すると、 \[ \frac{(n-1)s^2}{\chi^2_2} < \sigma^2 < \frac{(n-1)s^2}{\chi^2_1} \] という母分散 \(\sigma^2\) の \((1-\alpha) \times 100 \%\) 信頼区間が得られる。
なお、母平均 \(\mu\) が既知のときは、すでに述べたように、母分散は \[ \hat{\sigma^2} = \frac{1}{n}\sum_i(X_i - \mu)^2 \] で推定され、この標本分布は自由度 \(n\) の \(\chi^2\) 分布で記述されるので、母分散の信頼区間は、 \[ \frac{n \hat{\sigma^2}}{\chi^2_{n, 1-\alpha/2}} < \sigma^2 < \frac{n\hat{\sigma^2}}{\chi^2_{n, \alpha/2}} \] になる。なお、標準偏差の信頼区間は、分散の信頼区間の平方根で求められる。
ベニテングダケ(Amanita muscaria)は鮮やかな赤の傘を持つキノコで、傘に点在する白いいぼ、白いえら、同心円状の渦巻きなどの特徴がある。ベニテングダケの胞子は、白く、楕円形であり、7から 13μmの直径をもつ。 https://commons.wikimedia.org/wiki/File:Amanita_muscaria_3_vliegenzwammen_op_rij.jpg
いま、51のベニテングダケについて胞子の直径 \(X\) を計測したところ、以下のデータを得た。
# ベニテングダケ胞子の直径
dm <- c(10, 11, 12, 9, 10, 11, 13, 12, 10, 11, 11, 13, 9, 10, 9,
10, 8, 12, 10, 11, 9, 10, 7, 11, 8, 9, 11, 11, 10, 12, 10, 8, 7,
11, 12, 10, 9, 10, 11, 10, 8, 10, 10, 8, 9, 10, 13, 9, 12, 9, 9)
計測値は正規分布からの標本であると仮定したとき、以下の問に答えよ。
xbar <- mean(dm); xbar # 標本平均 10.09804
## [1] 10.09804
s <- sd(dm); s # 標本標準偏差 1.473159
## [1] 1.473159
n <- length(dm); n # 標本の大きさ 51
## [1] 51
# 1. 母平均の検定
mu0 <- 11
tv <- sqrt(n)*(xbar - mu0)/s; tv
## [1] -4.372434
tv <- abs(tv) # t 値
alpha1 <- 0.05; alpha2 <- 0.01
t5 <- qt(1 - alpha1/2, df = n - 1); t5 # 両側5%点
## [1] 2.008559
t1 <- qt(1 - alpha2/2, df = n - 1); t1 # 両側1%点
## [1] 2.677793
pv <- 2*(1 - pt(tv, df= n - 1)); pv # p 値
## [1] 6.224733e-05
\(t\)値は4.37で自由度50の\(t\)分布の両側1%点である2.678より大きいので1%有意であり、帰無仮説は棄却され、この標本の母平均 10.1μm は 11μm と有意に異なる。 また\(p\)値は\(\approx 0\)で非常に小さいので、帰無仮説は棄却される。なお、Rの関数 t.test()を用いても良い。
t.test(dm, mu = 11)
##
## One Sample t-test
##
## data: dm
## t = -4.3724, df = 50, p-value = 6.225e-05
## alternative hypothesis: true mean is not equal to 11
## 95 percent confidence interval:
## 9.683707 10.512372
## sample estimates:
## mean of x
## 10.09804
# 2. 母分散の検定
alpha <- 0.05 # 両側 5 %検定
sig02 <- 2
s^2 # 標本分散
## [1] 2.170196
x2 <- (n - 1)*s^2/sig02; x2 # カイ2乗値
## [1] 54.2549
chi1 <- qchisq(alpha/2, df = n - 1); chi1 # 下側 2.5%点
## [1] 32.35736
chi2 <- qchisq(1 - alpha/2, df = n - 1); chi2 # 上側 2.5%点
## [1] 71.4202
\(\chi^2\)値54.25で、これは自由度50の \(\chi^2\) 分布の下側2.5%点である32.36と上側2.5%点である71.42との間にあるので、5%有意でない。よって、帰無仮説は棄却されず母分散は2であると考えても良い。なお、1標本の分散に関する検定はパッケージEnvStatsのvarTest関数を用いて行うことができる。
EnvStats::varTest(dm, alternative = "two.sided", conf.level = 0.95, sigma.squared = 2)
##
## Results of Hypothesis Test
## --------------------------
##
## Null Hypothesis: variance = 2
##
## Alternative Hypothesis: True variance is not equal to 2
##
## Test Name: Chi-Squared Test on Variance
##
## Estimated Parameter(s): variance = 2.170196
##
## Data: dm
##
## Test Statistic: Chi-Squared = 54.2549
##
## Test Statistic Parameter: df = 50
##
## P-value: 0.6310443
##
## 95% Confidence Interval: LCL = 1.519315
## UCL = 3.353481
# 3. 母平均の95%信頼区間 (95% C.I.)
alpha <- 0.05
tquantile <- qt(1 - alpha/2, df = n - 1)
LB <- xbar - tquantile * s / sqrt(n)
UB <- xbar + tquantile * s / sqrt(n)
c(LB, UB) # 母平均μの95%信頼区間
## [1] 9.683707 10.512372
これより、母平均の95%信頼区間は、\(9.68 < \mu < 10.51\) となった。\(\mu_0 = 11\) はこの信頼区間に含まれていないので、5%有意であることが分かる。なお、t.test関数を用いた場合の結果にも母平均の95%信頼区間が示されている。
# 4. 母分散の90%信頼区間 (90% C.I.)
alpha <- 0.10
tchisq.1 <- qchisq(1 - alpha/2, df = n - 1)
tchisq.2 <- qchisq(alpha/2, df = n - 1)
LB <- (n - 1) * s^2 / tchisq.1
UB <- (n - 1) * s^2 / tchisq.2
c(LB, UB) # 母分散の90%信頼区間
## [1] 1.607438 3.121304
sqrt(c(LB, UB)) # 母標準偏差σの90%信頼区間
## [1] 1.267848 1.766721
これより、母分散と母標準偏差の90%信頼区間はそれぞれ、\(1.60 < \sigma^2 < 3.12、1.27 <\sigma < 1.77\) となることが分かる。なお、EnvStatsのvarTest関数を用いて行う場合は、conf.levelを0.90とする。
EnvStats::varTest(dm, alternative = "two.sided", conf.level = 0.90, sigma.squared = 2)
##
## Results of Hypothesis Test
## --------------------------
##
## Null Hypothesis: variance = 2
##
## Alternative Hypothesis: True variance is not equal to 2
##
## Test Name: Chi-Squared Test on Variance
##
## Estimated Parameter(s): variance = 2.170196
##
## Data: dm
##
## Test Statistic: Chi-Squared = 54.2549
##
## Test Statistic Parameter: df = 50
##
## P-value: 0.6310443
##
## 90% Confidence Interval: LCL = 1.607438
## UCL = 3.121304