このページに対応するRmdファイル：GitHub

1 Generating random numbers / 乱数の生成

確率分布 \(F\) に従う確率変数 \(X\) があったとき，この確率変数の実現値を疑似的に作り出す手続きを乱数の生成と呼ぶ．

\[ X_1, \ldots, X_n \sim F \]

疑似乱数は rxxx 関数で生成することができる．

  • xxx は確率分布の略称（一様分布 unif, 正規分布 norm, 二項分布 binom）

たとえば，一様分布から乱数をドローするには，runif関数を使う．この関数名は r (random) + unif (uniform distribution) から来ている．

runif(n = サンプルサイズ, min = 下限, max = 上限) で乱数を生成．

runif(n = 5, min = 0, max = 1)
## [1] 0.77570328 0.83649534 0.26750027 0.03175042 0.86424274
runif(5, 0, 1)  # 引数の順番を守れば「n=」などを省略可
## [1] 0.34332467 0.09578165 0.22716778 0.74271783 0.73830497
runif(5)  # min = 0 かつ max = 1 の場合は省略可（デフォルトの設定）
## [1] 0.3403366 0.6556750 0.5396952 0.2703124 0.8193960

乱数をドローするたびに異なる値が出力される．

1.1 Set seed

再現性のために乱数のシードを指定することができる．

set.seed(0)
runif(5)
## [1] 0.8966972 0.2655087 0.3721239 0.5728534 0.9082078
set.seed(0)
runif(5)
## [1] 0.8966972 0.2655087 0.3721239 0.5728534 0.9082078

2 Probability and Probability distribution / 確率と確率分布

2.1 Bernoulli distribution / ベルヌーイ分布

\[ \Pr(X = 1) = p, \quad \Pr(X = 0) = 1-p \]

例：60%の確率で表が出るコインをトスし，表が出れば1，裏が出れば0とする．

rbinom(n = 1, size = 1, prob = 0.6)
## [1] 1
table(rbinom(n = 10000, size = 1, prob = 0.6))
## 
##    0    1 
## 3969 6031

確率分布（確率変数が取りうる値とその確率を対応させたもの）は次の通り．

barplot(c(0.4, 0.6), names.arg = 0:1, ylab = "Probability")

2.2 Binomial distribution / 二項分布

\[ \Pr(X = k) = \left( \matrix{n \\ k} \right) p^k (1-p)^{n-k} \]

例：60%の確率で表が出るコインを5回トスしたときに表が出る回数．

rbinom(n = 1, size = 5, prob = 0.6)
## [1] 4
table(rbinom(n = 10000, size = 5, prob = 0.6))
## 
##    0    1    2    3    4    5 
##   90  790 2299 3468 2578  775

上の式より，表が 0, 1, 2, 3, 4, 5 回出る理論上の確率は次のように計算される．

binom_prob <- c(
  choose(5, 0) * 0.6^0 * (1 - 0.6)^(5 - 0),  # choose(n, k) = nCk
  choose(5, 1) * 0.6^1 * (1 - 0.6)^(5 - 1),
  choose(5, 2) * 0.6^2 * (1 - 0.6)^(5 - 2),
  choose(5, 3) * 0.6^3 * (1 - 0.6)^(5 - 3),
  choose(5, 4) * 0.6^4 * (1 - 0.6)^(5 - 4),
  choose(5, 5) * 0.6^5 * (1 - 0.6)^(5 - 5))
binom_prob
## [1] 0.01024 0.07680 0.23040 0.34560 0.25920 0.07776

よって，確率分布は次の通り．

barplot(binom_prob, names.arg = 0:5, ylab = "Probability")

2.3 Uniform distribution / 一様分布

\[ X \sim U(\mbox{min}, \mbox{max}) \]

runif(n = 1, min = 0, max = 1)
## [1] 0.117898
ru <- runif(n = 10000, min = 0, max = 1)
hist(ru, breaks = 100)

2.3.1 Prob

連続確率変数について所与の区間の値を取る確率は密度関数 \(f\) を用いて次で与えられる．

\[ \Pr(a \le X \le b) = \int_a^b f(x) dx = \int_a^b \frac{1}{\mbox{max} - \mbox{min}} dx \]

最大値1，最小値0の一様分布の場合は次のとおり（ただし \(0<a<b<1\)）．

\[ \Pr(a \le X \le b) = \int_a^b \frac{1}{1-0} dx = b-a \]

たとえば \(\Pr(0.2 \le X \le 0.9) = 0.7\) である．

head(data.frame(random_num = ru, is_in_range = ru >= 0.2 & ru <= 0.9))
##   random_num is_in_range
## 1  0.4772929        TRUE
## 2  0.7689994        TRUE
## 3  0.6008625        TRUE
## 4  0.2997955        TRUE
## 5  0.6545341        TRUE
## 6  0.3811610        TRUE
table(ru >= 0.2 & ru <= 0.9)
## 
## FALSE  TRUE 
##  2944  7056

分布関数 \(F\) に基づいて下図のように考えてもよい．

この場合，所与の \(x\) の値に対応する累積確率を返す関数 punif を用いて次のように計算することもできる． 関数名は p (probability) + unif より．

punif(0.9) - punif(0.2)
## [1] 0.7

つまり，punif 関数と密度関数 \(f\) および分布関数 \(F\) は次のような関係にある．

\[ \mbox{punif}(x, \mbox{min}, \mbox{max}) = \int_{\mbox{min}}^{x} f(u) du = F(x) = \Pr(X \le x) \] \[ \mbox{where} \quad f(u) = \frac{1}{\mbox{max} - \mbox{min}} \]

2.4 Normal distribution / 正規分布

\[ X \sim N(\mu, \sigma^2) \]

rnorm(n = 1, mean = 0, sd = 1)
## [1] 0.9210962
hist(rnorm(n = 10000, mean = 00, sd = 1), breaks = 100)

2.4.1 Prob

平均 0 分散 1 の正規分布を標準正規分布と呼ぶ．

\[ Z \sim N(0, 1) \]

標準正規分布の密度関数 \(\phi\) について \(z = 1.96\) より左側の面積は確率 \(\Pr(Z \le 1.96)\) である．

\[ \Pr(Z \le 1.96) = \int_{-\infty}^{1.96} \phi (x) dx \]

一様分布のときと同様に，累積分布関数に基づいて次のように計算できる．

pnorm は与えられた \(z\) スコアに対応する累積確率を返す関数． 関数名は p (probability) + norm (normal distribution) から来ている．

pnorm(q = 1.96, mean = 0, sd = 1)
## [1] 0.9750021

次のように正規乱数を用いて確認することもできる．

rn <- rnorm(n = 10000, mean = 0, sd = 1)
head(data.frame(random_num = rn, is_in_range = rn <= 1.96))
##   random_num is_in_range
## 1  0.9134180        TRUE
## 2 -2.8988797        TRUE
## 3 -0.0213748        TRUE
## 4 -0.9831606        TRUE
## 5  0.1166303        TRUE
## 6 -0.1551706        TRUE
table(rn <= 1.96)
## 
## FALSE  TRUE 
##   265  9735

2.4.2 qnorm

pnorm とは逆に，qnorm は与えられた確率に対応する \(z\) スコアを返す関数． 

qnorm(p = 0.975, mean = 0, sd = 1)
## [1] 1.959964
qnorm(p = 0.975)  # for N(0,1), the mean and sd arguments may be omitted 
## [1] 1.959964

2.4.3 dnorm

確率密度． dnorm = d (density) + norm (normal distribution).

dnorm(x = 1.96, mean = 0, sd = 1)
## [1] 0.05844094

標準正規分布の確率密度 \(\phi\) を積分する（確率密度関数とX軸の間の面積を計算する）と，積分区間の間の値を取る確率が得られる．

\[ \Pr(-1.96 \le Z \le 1.96) = \int_{-1.96}^{1.96} \phi(z) dz \approx 0.95 \]

integrate(dnorm, -1.96, 1.96)  # 積分: 約 95%
## 0.9500042 with absolute error < 1e-11

3 Expectation and Variance of random variable

離散確率変数 \(X\) の期待値と分散は確率関数 \(f\) を用いて以下のように定義される．

\[ \mu = E(X) = \sum_j^k x_j f(x_j) \]

\[ \sigma^2 = V(X) = E[(X-E(X))^2] = \sum_j^k (x_j - E(X))^2 f(x_j) = E(X^2) - (E(X))^2 \]

連続確率変数の場合（ここでは \(f\) は確率密度関数）

\[ \mu = E(X) = \int_{-\infty}^{\infty} x f(x) dx \]

\[ \sigma^2 = V(X) = \int_{-\infty}^{\infty} (x - \mu)^2 f(x) dx\]

3.1 Dice roll

例：サイコロの出目 \(X\) の期待値と分散 / Example of dice roll

\[ E(X) = \sum_{j=1}^6 x_j \frac{1}{6} = 1 \times \frac{1}{6} + 2 \times \frac{1}{6} + 3 \times \frac{1}{6} + 4 \times \frac{1}{6} + 5 \times \frac{1}{6} + 6 \times \frac{1}{6} = \frac{7}{2} \]

dice <- 1:6  # c(1, 2, 3, 4, 5, 6) と同じ
sum(dice * (1/6))  # definition of expectation
## [1] 3.5
mean(dice)  # mean
## [1] 3.5

\[ V(X) = \sum_{j=1}^6 (x_j - E(X))^2 \frac{1}{6} = \left( 1-\frac{7}{2} \right) ^2 \frac{1}{6} + \left( 2-\frac{7}{2} \right) ^2 \frac{1}{6} + \cdots + \left( 6-\frac{7}{2} \right) ^2 \frac{1}{6} \approx 2.917 \]

sum((dice - mean(dice))^2 * (1/6))  # variance: sum[(x-E(X))*f(x)]
## [1] 2.916667

なお，var(dice) は 3.5 となり，上記と一致しないが，これは var 関数が標本を引数として標本分散（不偏分散）を計算するためである． dice は母集団から抽出された標本ではないため，標本分散として確率変数の分散を求めることはできない．

\(V(X) = E(X^2) - (E(X))^2\) でもある．

sum(dice^2 * (1/6)) - sum(dice * (1/6))^2  # variance: E(X^2)-E(X)^2
## [1] 2.916667

3.2 Uniformly distributed random variable

一様分布 \(U(a, b)\) から乱数を生成して標本平均と標本分散を計算し，それが理論値と一致するかを確認する．

\[ E(X) = \int_a^b x \frac{1}{b-a} dx = \frac{1}{b-a} \left[ \frac{x^2}{2} \right]_a^b = \frac{1}{b-a} \times \frac{b^2 - a^2}{2} = \frac{a+b}{2} \]

\[ V(X) = \int_a^b (x - \mu)^2 \frac{1}{b-a} dx = \frac{1}{b-a} \left[ \frac{ \left( x - \frac{a+b}{2} \right)^3}{3} \right]_a^b = \frac{(b - a)^2}{12} \]

\(X \sim U(-1, 1)\) とすると，期待値と分散は次の値となる．

\[ E(X) = \frac{a+b}{2} = \frac{-1+1}{2} = 0 \]

\[ V(X) = \frac{(b - a)^2}{12} = \frac{(1 - (-1))^2}{12} = \frac{1}{3} \]

「サンプルサイズ 100 で一様乱数を生成して平均と分散を計算する」を 1,000 回繰り返し，1,000個の平均と分散の分布を観察する． このようなシミュレーションを「モンテカルロ・シミュレーション」と呼ぶ．

  • for 文を使って，乱数生成，平均と分散の計算，記録を行う
    • 計算結果を格納するオブジェクトを runif_mean と runif_var という名前で作成．NULL を代入することで，空のオブジェクトを作成
    • c(直前までの計算結果のベクトル, 今回の計算結果) で append
runif_mean <- NULL  # 標本平均を格納しておくオブジェクト
runif_var <- NULL  # 標本分散を 〃
for (i in 1:1000) {
  set.seed(i)  # 再現性のため
  # 乱数を生成
  runif_i <- runif(n = 100, min = -1, max = 1)
  # 平均・分散を計算して上で作成したオブジェクトに格納 (append) 
  runif_mean <- c(runif_mean, mean(runif_i))
  runif_var <- c(runif_var, var(runif_i))
}
summary(runif_mean)
##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## -1.842e-01 -4.050e-02 -4.024e-05 -1.030e-03  3.643e-02  2.225e-01
summary(runif_var)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.2413  0.3127  0.3327  0.3333  0.3545  0.4225

4 Law of large numbers / 大数の法則

ラフに言えば，サンプルサイズが大きくなれば標本平均が母平均に収束するという法則．

確率変数 \(X_1, \ldots, X_n\) が互いに独立に同一の分布に従い \(E(X_i) = \mu < \infty\) のとき，任意の \(\epsilon > 0\) に対して次が成立（証明はチェビシェフの不等式によるものが有名）．

\[ \Pr ( |\bar{X} - \mu| \le \epsilon ) \to 1 \quad (\mbox{as } n \to \infty) \]

正規分布に従う確率変数 \(X_i \sim N(0, 1)\) の標本平均 \(\bar{X}_n = (1/n) \sum_{i=1}^n X_n\) がサンプルサイズ \(n\) とともにどのように変化するかを観察する．

rn <- rn_mean <- NULL
for (i in 1:10000) {
  rn <- c(rn, rnorm(n = 1)) 
  rn_mean <- c(rn_mean, mean(rn)) 
}
head(cbind(random_num = rn, cumulative_mean = rn_mean))
##      random_num cumulative_mean
## [1,]  0.2517114      0.25171138
## [2,] -1.0946937     -0.42149116
## [3,]  0.3976428     -0.14844649
## [4,] -0.9963020     -0.36041037
## [5,]  0.1005780     -0.26821269
## [6,]  0.9536803     -0.06456386
plot(rn_mean, type = "l", log = "x")
abline(h = 0, col = "grey", lty = 1)  # mu = 0 (true mean) 
abline(h = c(0.05, -0.05), col = "grey", lty = 2)  # epsilon = 0.05

余談だが，実は次のような短いコードでも同じ（出力は省略）．

plot(cumsum(rnorm(n = 10000)) / 1:10000, type = "l", log = "x")

5 Central limit theorem / 中心極限定理

平均 \(\mu\) 分散 \(\sigma^2\) の母集団から得られた標本の平均 \(\bar{X}\) は次の正規分布に従う．

\[ \bar{X} \sim N \left( \mu, \frac{\sigma^2}{n} \right) \]

ポイント

  • 一つの標本から計算される標本平均は特定の値 \(\bar{x}\) をとるが，同じ母集団から何度もサンプリングを行いその度に繰り返し標本平均を計算すればその標本平均はばらつく．これを確率変数 \(\bar{X}\) として扱う．
  • サンプルサイズ \(n\) が大きくなると標本平均は母平均に近い値を取る．
  • 母集団分布によらずに（例外あり）標本平均は正規分布に従う．
  • サンプルサイズ \(n\) がある程度大きくないと正規分布に近似されない場合がある．
  • 中心極限定理は母平均の信頼区間の計算や平均値の差の検定を行う上で重要．

以下では一様確率変数の標本平均の分布を求める．標本サイズが大きくなるほど正規分布に近づき，分布の分散が小さくなることが分かる．

5.1 \(n=1\)

これは乱数そのものの分布と等しい．

\[ \bar{X} = X_i, \quad X_i \sim U(0, 1) \]

x_bar <- NULL  # create empty (null) object
for (i in 1:10000) {
  x_bar_i <- mean(runif(n = 1, min = 0, max = 1))
  x_bar <- c(x_bar, x_bar_i)
}
hist(x_bar, breaks = 50, main = "Sample mean of uniform random variables with sample size n=1")

var(x_bar)  # シミュレーションによる分散
## [1] 0.08302284
(1/12) / 1  # 理論上の分散
## [1] 0.08333333

5.2 \(n=2\)

\[ \bar{X} = \frac{1}{2} \sum_{i=1}^{2} X_i, \quad X_i \sim U(0, 1) \]

x_bar <- NULL
for (i in 1:10000) {
  x_bar_i <- mean(runif(n = 2, min = 0, max = 1))
  x_bar <- c(x_bar, x_bar_i)
}
hist(x_bar, breaks = 50, main = "Sample mean of uniform random variables with sample size n=2")

var(x_bar)  # シミュレーションによる分散
## [1] 0.04145714
(1/12) / 2  # 理論上の分散
## [1] 0.04166667

5.3 \(n=10\)

\[ \bar{X} = \frac{1}{10} \sum_{i=1}^{10} X_i, \quad X_i \sim U(0, 1) \]

x_bar <- NULL
for (i in 1:10000) {
  x_bar_i <- mean(runif(n = 10, min = 0, max = 1))
  x_bar <- c(x_bar, x_bar_i)
}
hist(x_bar, breaks = 50, main = "Sample mean of uniform random variables with sample size n=10",
     xlim = c(0, 1))

var(x_bar)  # シミュレーションによる分散
## [1] 0.008489861
(1/12) / 10  # 理論上の分散
## [1] 0.008333333

5.4 For advanced students / 中級者向け

for 文を使わずに sapply 関数を使って次のように書くこともできる： 

hist(sapply(X = 1:10000, FUN = function (x) mean(runif(10)) ), breaks = 50)

（出力は省略）

5.5 Normal distribution, \(n=2\)

正規母集団の場合はサンプルサイズが小さくても正規分布に従う（再生性と呼ばれる性質）．

\[ \bar{X} = \frac{1}{2} \sum_{i=1}^2 X_i, \quad X_i \sim N(0, 1) \]

x_bar <- NULL
for (i in 1:10000) {
  x_bar_i <- mean(rnorm(n = 2, mean = 0, sd = 1))
  x_bar <- c(x_bar, x_bar_i)
}
hist(x_bar, breaks = 50)

var(x_bar)  # シミュレーションによる分散
## [1] 0.4913276
1 / 2  # 理論上の分散
## [1] 0.5

この性質は，t検定による平均値の差の検定を適用してよいか判断する際の次の条件と関連する．

  • 母集団が正規分布であるか分からないがサンプルサイズが大きい → 中心極限定理により \(\bar{X} \sim N\)
  • 母集団が正規分布 → 正規分布の再生性により \(\bar{X} \sim N\)

5.6 Log-normal distribution, \(n=10\)

母集団分布の形状によってはサンプルサイズが小さいと正規分布で十分に近似できない． たとえば，対数正規分布のような歪んだ分布からサンプリングされた標本の場合は次のようになる．

\[ \bar{X} = \frac{1}{10} \sum_{i=1}^{10} X_i, \quad X_i \sim \ln N(0, 1) \]

x_bar <- NULL
for (i in 1:10000) {
  x_bar_i <- mean(rlnorm(n = 10, meanlog = 0, sdlog = 1))
  x_bar <- c(x_bar, x_bar_i)
}
hist(x_bar, breaks = 50)

var(x_bar)  # シミュレーションによる分散
## [1] 0.4735483
(exp(2) - exp(1)) / 10  # 理論上の分散
## [1] 0.4670774

6 Confidence interval / 信頼区間

パラメタ \(\theta\) に関する信頼水準 \((1-\alpha)\) の信頼区間 \([\hat{\theta}_\mbox{lower}, \hat{\theta}_\mbox{upper}]\) は，直感的には次のように理解できる．

\[ \Pr(\hat{\theta}_\mbox{lower} \le \theta \le \hat{\theta}_\mbox{upper}) = 1 - \alpha \]

たとえば，母平均 \(\mu\) の95%信頼区間は \(\sigma\) が既知の場合次のように求められる．

\[ \Pr \left( -1.96 \le \frac{\bar{X} - \mu}{\sigma/\sqrt{n}} \le 1.96 \right) = 0.95 \] \[ \therefore \Pr \left( \bar{X} - 1.96 \times \frac{\sigma}{\sqrt{n}} \le \mu \le \bar{X} + 1.96 \times \frac{\sigma}{\sqrt{n}} \right) = 0.95 \]

6.1 Population mean with known \(\sigma\) / 母集団の標準偏差が既知の場合の母平均

標本平均が100，「母集団」の標準偏差が2，サンプルサイズが100の場合，信頼水準95%の信頼区間は次のように計算される．

\[ \bar{x} \pm z \frac{\sigma}{\sqrt{n}} = 100 \pm 1.96 \frac{2}{\sqrt{100}} = 100 \pm 0.392, \quad \mbox{95% CI}: [99.61, 100.39] \]

\(z = 1.96\) は，標準正規分布に従う確率変数 \(Z\) について \(\Pr(Z \le z) = 0.975\) を満たす値．

# qnorm(p = 0.025, mean = 0, sd = 1)  # P(Z<z) = 0.025
qnorm(p = 0.975, mean = 0, sd = 1)  # P(Z<z) = 0.975
## [1] 1.959964

100 - qnorm(p = 0.975) * 2 / sqrt(100)  # lower limit
## [1] 99.60801
100 + qnorm(p = 0.975) * 2 / sqrt(100)  # upper limit
## [1] 100.392

6.1.1 Simulation

信頼区間とは，同じ母集団から何度も繰り返し「サンプリング → 母平均の信頼区間を計算」をしたときに，母平均が \((1-\alpha)\) の割合で含まれる区間である．

このことを，次のようなシミュレーションで確認してみよう．

  • 母集団 … 母平均 \(\mu = 100\), 母分散 \(\sigma^2 = 2^2=4\) の正規分布
  • 標本と信頼水準 … サンプルサイズ \(n = 100\), \(\alpha = 0.05\)
ci_list <- NULL
for (i in 1:1000) {
  sample_i <- rnorm(n = 100, mean = 100, sd = 2)
  ci_i_lower <- mean(sample_i) - (1.96 * 2 / sqrt(100))
  ci_i_upper <- mean(sample_i) + (1.96 * 2 / sqrt(100))
  ci_list <- rbind(ci_list, c(lower = ci_i_lower, upper = ci_i_upper))
}
ci_list <- data.frame(ci_list, is_incl_mu = ci_list[,1] <= 100 & ci_list[,2] >= 100)
head(ci_list)
##      lower     upper is_incl_mu
## 1 99.17879  99.96279      FALSE
## 2 99.72957 100.51357       TRUE
## 3 99.45163 100.23563       TRUE
## 4 99.88772 100.67172       TRUE
## 5 99.40413 100.18813       TRUE
## 6 99.60410 100.38810       TRUE
table(ci_list$is_incl_mu)
## 
## FALSE  TRUE 
##    43   957

6.2 Population mean with unknown \(\sigma\) / 母集団の標準偏差が未知の場合の母平均

標本平均が100，「標本」標準偏差が2，サンプルサイズが100の場合，信頼水準95%の信頼区間は次のようになる．

\[ \bar{x} \pm t \frac{s}{\sqrt{n}} = 100 \pm 1.98 \frac{2}{\sqrt{100}} = 100 \pm 0.396, \quad \mbox{95% CI}: [99.60, 100.40] \]

qt(p = 0.975, df = 100-1)  # P(T<t) = 0.025
## [1] 1.984217
100 - qt(p = 0.975, df = 100-1) * 2 / sqrt(100)  # lower limit
## [1] 99.60316
100 + qt(p = 0.975, df = 100-1) * 2 / sqrt(100)  # upper limit
## [1] 100.3968

6.2.1 Using t.test function

\(t\)-検定を行うt.testという関数は信頼区間をついでに出力してくれる．

まずは t.test 関数を使わない場合． 信頼水準95%とする．

wage <- c(1000, 1200, 1300, 1200, 1150, 1000, 1450, 1500, 1150, 1350)
c(mean(wage), sd(wage), length(wage))  # mean, standard deviation, sample size
## [1] 1230.0000  170.2939   10.0000
mean(wage) - qt(p = 0.975, df = length(wage)-1) * sd(wage) / sqrt(length(wage))  # lower bound 
## [1] 1108.179
mean(wage) + qt(p = 0.975, df = length(wage)-1) * sd(wage) / sqrt(length(wage))  # upper bound 
## [1] 1351.821

同じ結果が t.test 関数で確認できる．

t.test(wage)  # default: 95% CI
## 
##  One Sample t-test
## 
## data:  wage
## t = 22.841, df = 9, p-value = 2.806e-09
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  1108.179 1351.821
## sample estimates:
## mean of x 
##      1230

信頼水準 (\(1 - \alpha\)) はデフォルトで95%に設定されている． conf.level 引数を指定すればそれ以外の信頼水準に変更できる．

t.test(wage, conf.level = 0.99)  # 99% CI
## 
##  One Sample t-test
## 
## data:  wage
## t = 22.841, df = 9, p-value = 2.806e-09
## alternative hypothesis: true mean is not equal to 0
## 99 percent confidence interval:
##  1054.991 1405.009
## sample estimates:
## mean of x 
##      1230

6.3 Population proportion / 母比率

標本の比率 \(p\)，サンプルサイズ \(n\)，z 値 (such that \(\Pr(Z \le z) = 1 - \alpha\)) より次のように母比率の信頼区間が計算される．

\[ p \pm z \sqrt{\frac{p (1-p)}{n}} \]

Example of Banerjee, Duflo, and Glennerster (BMJ 2010)

382人の子供のうち148人がワクチンの予防接種を受けた． 

n1 <- 148; n0 <- 382  # intervention B
p <- n1 / n0
p
## [1] 0.3874346
p - qnorm(p = .975) * sqrt(p * (1 - p) / n0)  # lower limit of 95% CI
## [1] 0.3385815
p + qnorm(p = .975) * sqrt(p * (1 - p) / n0)  # upper limit
## [1] 0.4362876

注：この数値は論文中で報告されている数値とは異なる． 地理的なブロック（村）を用いてクラスター無作為化対照試験が実行されたため，著者らはこのデータの階層性を考慮したより複雑な計算をしていると思われる．

7 Statistical hypothesis testing / 統計的仮説検定

仮説検定の手続き：

  • 帰無仮説 \(H_0\) と対立仮説 \(H_1\) を宣言
    • 帰無仮説は母集団パラメタが特定の値（または他のパラメタ）と等しいという形で定式化され［例：\(\theta = \theta_0\)］，対立仮説は特定の値（または他のパラメタ）と異なるという形で定式化される［例：\(\theta > \theta_0\) または \(\theta \ne \theta_0\)］
  • 検定統計量を選択
    • 例：平均値の差の検定で母分散が既知であれば \(z\)，平均値の差の検定で母分散が未知であれば \(t\)
  • 帰無仮説のもとでの検定統計量の標本分布を導出
    • 例：検定統計量が \(z\) であれば標準正規分布，検定統計量が \(t\) であれば \(t\) 分布
  • 有意水準 \(\alpha\) を選択
    • \(\alpha = 0.05\) または \(\alpha = 0.01\) が慣習的によく用いられる
  • \(p\) 値を計算
    • \(p\) 値とは，検定統計量が実際に観察された値あるいはそれ以上に極端な値を取る確率として定義される
    • 対立仮説が \(\theta > \theta_0\) の形の場合は片側 \(p\) 値を計算
    • 対立仮説が \(\theta \ne \theta_0\) の形の場合は両側 \(p\) 値を計算
  • \(p \le \alpha\) ならば帰無仮説を棄却し対立仮説を採択する
    • 注：\(p > \alpha\) のときに「帰無仮説を採択する」ことはできない

7.1 One-sample \(t\) test for population mean / 1標本の母平均の \(t\) 検定

帰無仮説：

\[ H_0: \mu = \mu_0 \]

たとえば，標本平均 \(\bar{x} = 1230\) に対して帰無仮説が \(H_0: \mu = 1100\) であった（\(\bar{x} > \mu_0\)）とすると，我々が計算したい \(p\) 値は次のように定義できるはず（あくまでイメージ）．

\[ p = \Pr(\bar{X} \ge 1230) \quad \mbox{or} \quad p = 2 \times \Pr(\bar{X} \ge 1230) \]

つまり，この確率を計算するためには \(\bar{X}\) すなわち標本平均という確率変数が従う分布を知る必要がある．

  • 母集団分布は未知だがサンプルサイズが十分に大きい → 中心極限定理より \(\bar{X} \sim N\)
    • 母分散が既知であれば \(N\)，未知であれば \(t\) が検定統計量が従う分布

\[ Z = \frac{\bar{X} - \mu_0}{\sigma/\sqrt{n}} \sim N(0, 1), \quad T = \frac{\bar{X} - \mu_0}{s/\sqrt{n}} \sim t(n-1) \]

  • 母集団が正規分布 → 正規分布の再生性より \(\bar{X} \sim N\)
    • 同上

\(N(0,1)\) や \(t(n-1)\) は，帰無仮説が正しいもとで検定統計量が従う分布である．

たとえば，帰無仮説 \(H_0: \mu = 1100\) が実際に正しければ，標本平均は 1100 に近い値をとる可能性が高く，標本平均の期待値は 1100 になるはずである． そこで，\(N(1100,1)\) から乱数を生成してサイズ 5 の標本を1万個作り，その標本から得られる検定統計量がどのような分布を形成するかを観察してみよう． 母分散 \(\sigma^2\) は未知であるとすると，理論上は \(t\) 分布が出現する．

n <- 5
t_value <- NULL
for (i in 1:10000) {
  rn <- rnorm(n = n, mean = 1100)
  t_value <- c(t_value, (mean(rn)-1100) / (sd(rn)/sqrt(n)))
}
hist(t_value, breaks = 200, probability = TRUE, xlim = c(-6, 6), ylim = c(0, 0.4))
# 理論上の分布 t(n-1) を重ね描き
x_vals <- seq(-4, 4, length.out = 100)
lines(x_vals, dt(x_vals, df = n-1), col = "blue", lwd = 2)  # t(n-1)
lines(x_vals, dnorm(x_vals), col = "red", lwd = 1, lty = 2)  # N (for reference)
legend(x = "topleft", legend = c("t(n-1)", "N(0,1)"), col = c("blue", "red"), lty = 1:2, cex = 1.5)

もし実際に観察された検定統計量の値がこの分布の端（たとえば \(t=2.4\)）に位置するのであれば，それは「帰無仮説が正しいときに得られるはずの値から大きく乖離した平均値が観察されていること」を意味する． このとき帰無仮説とデータは矛盾していると判断し，帰無仮説を棄却する．

帰無仮説とデータの乖離を測る指標として次の \(p\) 値が計算される．

\[ p = \Pr(T \ge 2.4) \quad \mbox{or} \quad p = 2 \times \Pr(T \ge 2.4) \]

累積分布関数を用いて汎用に書けば次のようになる．

\[ p = 1 - \Pr(T \le |t|) \quad \mbox{or} \quad p = 2 \times (1 - \Pr(T \le |t|)) \]

7.1.1 Calculation from the sample mean and sd of the data / データの標本平均・標準偏差から計算

賃金の母集団が正規分布に従っており，\(\sigma\) は未知，\(\bar{x} = 1230, s=170.3, n = 10, \alpha=0.05\) とする．

\[ H_0: \mu = 1100, \quad H_1: \mu \ne 1100 \]

\[ t = \frac{\bar{x} - \mu_0}{s/\sqrt{n}} = \frac{1230 - 1100}{170.3 / \sqrt{10}} \approx 2.41 \]

\[ \Pr (T \le 2.41) = 0.98 , \quad p = 2 \times (1 - \Pr (T \le 2.41)) \approx 0.039 \]

t_value <- (1230 - 1100) / (170.3 / sqrt(10))  # t value
t_value
## [1] 2.413952
pt(q = t_value, df = 10 - 1)  # Pr(T <= t)
## [1] 0.9805023
1 - pt(q = t_value, df = 10 - 1)  # one sided p value ... p = Pr(T > t)
## [1] 0.01949773
2 * (1 - pt(q = t_value, df = 10 - 1))  # two sided p value
## [1] 0.03899546

よって，\(H_0\) は棄却され \(H_1\) が採択される．

7.1.1.1 Supplement

直感的には（＝厳密ではないけれども）次のように理解できる．

我々は，データと帰無仮説がどの程度矛盾しているかを知りたい． もし矛盾の程度が大きければ，帰無仮説は間違っているだろうと考える（これは背理法のアイデア）． すなわち，この矛盾の程度の測定が仮説検定のポイントとなる．

矛盾の程度は \(p\) 値として測られる．これを測定するために次のように考える．

  • まず，注目するパラメタに対応する推定量は何で，その推定量（データの関数）がどのように分布するかを考える．たとえば，注目するパラメタが母平均 \(\mu\) であれば，それに対応する推定量は標本平均 \(\bar{X}\) である．サンプルサイズが大きければ中心極限定理より標本平均は正規分布に従う．
  • 次いで，帰無仮説が正しいときに，そのような母集団から得られる標本から計算される推定量はどのような分布に従うかを考える．たとえば，帰無仮説が \(\mu = 0\) でこれが実際に正しければ，その母集団から得られる標本の平均は 0 に近い値になるはずである．標本平均のばらつき（標準誤差）は母集団分布の分散およびサンプルサイズに依存するので，これらの情報を用いて標本平均をスケーリングする．これが検定統計量 \(Z = (\bar{X}-0)/(\sigma/\sqrt{n})\) となる．
  • この検定統計量は，標準正規分布に従う．これは「もし帰無仮説が本当に正しければ，検定統計量は標準正規分布に従うはずである」という意味であり，ラフに言えば「帰無仮説が正しければ，検定統計量は 0.5 や -1.0 のような 0 に近い値をとるはず．2.0 や -3.0 のように 0 から離れた値をとることはほとんどないはず．」ということである．
  • 以上に基づいて，「帰無仮説のもとで検定統計量が従う確率分布」において実際に観察された検定統計量がどこに位置するかを測る．このために計算されるのが \(p\) 値である（\(p = 2 \times (1 - \Pr(Z \le |z|))\)）．
    • 実際に観察された検定統計量が分布の中央に位置するなら \(p\) 値は大きな値を取り（例：\(p = 2 \times (1 - 0.5)=1\)），「帰無仮説と観察されたデータの間に矛盾は認められない」と結論付けられる．この場合，帰無仮説を棄却できない．言い換えれば「母平均が 0 と異なると結論付けることはできない」が結論となる．
    • 一方で実際に観察された検定統計量が分布の端に位置するなら \(p\) 値は小さな値を取り（例：\(p = 2 \times (1 - 0.99)=0.02\)），「帰無仮説と観察されたデータは矛盾する」と結論付けられる．この場合，帰無仮説は棄却され，対立仮説が採択される．言い換えれば「母平均は 0 と異なる値を取る」と結論付けることができる．ただしこの結論は母集団からのサンプリングという不確実性を含む値に基づいているため，ある一定の確率で誤る．すなわち「本当は帰無仮説が正しいのだが，偶然にも偏った値ばかりがサンプリングされてしまった結果，誤って帰無仮説が棄却された」（第一種の過誤）という可能性がある．この誤りの程度として許容できる水準は有意水準 \(\alpha\) として分析手続きにおいて明示的・定量的に扱われている．
  • いずれにせよ，上記の結論は我々分析者が想定する統計モデル（標本が同一母集団からの独立なサンプリングによって得られること，中心極限定理が適用できるような母集団でありサンプルサイズが十分に大きいこと or 正規母集団であること，母分散の値が正確であること，etc.）が正しいという前提のもとで正当化される．前提が間違っていれば結論は（実務的にはともかく，理論上は）意味を持たない．

7.1.2 Calculation from the data itself / データそのものから計算

上と同様に帰無仮説を \(\mu=1100\) とする．

wage <- c(1000, 1200, 1300, 1200, 1150, 1000, 1450, 1500, 1150, 1350)
c(mean(wage), sd(wage), length(wage))
## [1] 1230.0000  170.2939   10.0000
t_value <- (mean(wage) - 1100) / (sd(wage) / sqrt(length(wage)))
2 * (1 - pt(q = t_value, df = length(wage) - 1))  # two-sided p value
## [1] 0.0389899

t.test(ベクトル, mu = 帰無仮説の仮説値) 関数を使用．

t.test(wage, mu = 1100)  # default: two sided ... H1 is [mu != 1100]
## 
##  One Sample t-test
## 
## data:  wage
## t = 2.414, df = 9, p-value = 0.03899
## alternative hypothesis: true mean is not equal to 1100
## 95 percent confidence interval:
##  1108.179 1351.821
## sample estimates:
## mean of x 
##      1230

帰無仮説 \(\mu = \mu_0\) の \(\mu_0\) の値が異なったら \(p\) 値はどのように異なるか？

t.test(wage, mu = 100)$p.value
## [1] 5.942946e-09
t.test(wage, mu = 1000)$p.value
## [1] 0.002077277
t.test(wage, mu = 1100)$p.value
## [1] 0.0389899
t.test(wage, mu = 1200)$p.value
## [1] 0.5910512
t.test(wage, mu = 1230)$p.value
## [1] 1
t.test(wage, mu = 1250)$p.value
## [1] 0.718939

片側検定をする場合は alternative = "greater" または = "less" 引数で指定．

t.test(wage, mu = 1100, alternative = "greater")  # one sided ... H1 is [mu > 1100]
## 
##  One Sample t-test
## 
## data:  wage
## t = 2.414, df = 9, p-value = 0.01949
## alternative hypothesis: true mean is greater than 1100
## 95 percent confidence interval:
##  1131.284      Inf
## sample estimates:
## mean of x 
##      1230

このように信頼区間も同時に計算してくれる．

7.2 \(p\) and \(\alpha\)

\(p\) 値はどのように分布し，その値と有意水準 \(\alpha\) はどのような関係にあるのか？

「標準正規分布 \(N(0, 1)\) から \(n=100\) の乱数をドローして，その平均が 0 に等しいかどうかを検定する（\(H_0: \mu = 0\) とする一標本の母平均の検定）」という手続きを1,000回繰り返し，検定で計算されたp値の分布を観察する．

設定上，帰無仮説は棄却されないはずだが，現実には偶然に大きな値（または小さな値）ばかりがドローされることで標本平均が 0 から大きく乖離し，帰無仮説が棄却される場合がありうる．

t.test 関数で計算される p 値は t.test(...)$p.value で抽出できる．

t.test(x = 1:10)$p.value
## [1] 0.000278196
p_value <- NULL
for (i in 1:1000) {
  p_value_i <- t.test(rnorm(100, mean = 0, sd = 1))$p.value
  p_value <- c(p_value, p_value_i)
}
hist(p_value, breaks = 20, col = c("lightblue", rep("grey", 19)))

table(p_value <= 0.05)
## 
## FALSE  TRUE 
##   959    41

母集団分布の母平均は \(\mu = 0\) であるからこの帰無仮説は棄却されないことが望まれるが（\(p > \alpha\)），上記の通り5%程度の確率で誤って棄却されてしまう（\(p \le \alpha\), 第一種の過誤）．

7.2.1 \(t\)-test on a sample mean that is not normally distributed

中心極限定理が適用できないような標本平均に対して \(t\) 検定を適用したらどうなるだろうか？

歪んだ分布の例として対数正規分布 \(X \sim \ln N (\mu, \sigma^2)\) を取り上げ，ここからサイズ 10 の標本を作成する． 母平均が理論上の平均値 \(\exp(\mu + \sigma^2/2)\) と一致するかを検定する．

p_value <- numeric(1000)  # 事前にベクトルを作成すると計算速度が向上
for (i in 1:1000) {
  p_value[i] <- t.test(rlnorm(n = 10), mu = exp(0 + 1/2))$p.value
}
hist(p_value, breaks = 20, col = c("lightblue", rep("grey", 19)))

table(p_value <= 0.05)
## 
## FALSE  TRUE 
##   833   167

\(p \le \alpha\) （第一種の過誤）となる割合が \(\alpha\) より遥かに大きいことが確認できる．

7.3 Two-sample \(t\) test for population mean / 2標本の母平均の \(t\) 検定

7.3.1 Independent \(t\) test

\[ T = \frac{\bar{X}_1 - \bar{X}_2}{\sqrt{s_1^2 / n_1 + s_2^2 / n_2}} \sim t \left( \frac{[s_1^2 / n_1 + s_2^2 / n_2]^2}{\frac{(s_1^2 / n_1)^2}{n_1 - 1} + \frac{(s_2^2 / n_2)^2}{n_2 - 1}} \right) \]

7.3.1.1 Example of “wage” data

データ

wage_jp <- c(1000, 1200, 1300, 1200, 1150, 1000, 1450, 1500, 1150, 1350)  # Japan
wage_us <- c(900, 1300, 1200, 800, 1600, 850, 1000, 950)  # US
c(length(wage_jp), mean(wage_jp), sd(wage_jp))
## [1]   10.0000 1230.0000  170.2939
c(length(wage_us), mean(wage_us), sd(wage_us))
## [1]    8.0000 1075.0000  272.5541

定義に基づいて計算する場合．

t_value <- (mean(wage_jp) - mean(wage_us)) / 
  sqrt(var(wage_jp) / length(wage_jp) + var(wage_us) / length(wage_us))
df_wage <- (var(wage_jp) / length(wage_jp) + var(wage_us) / length(wage_us))^2 / 
  (var(wage_jp)^2 / length(wage_jp)^2 / (length(wage_jp)-1) +
     var(wage_us)^2 / length(wage_us)^2 / (length(wage_us) - 1))
c(t_value, df_wage, pt(q = t_value, df = df_wage))
## [1]  1.4041264 11.2050309  0.9063064
2 * (1 - pt(q = abs(t_value), df = df_wage))
## [1] 0.1873871

t.test(x = 1つ目のベクトル, y = 2つ目のベクトル) のように指定する．

t.test(wage_jp, wage_us)  # default: Welch's t test (assuming unequal variance)
## 
##  Welch Two Sample t-test
## 
## data:  wage_jp and wage_us
## t = 1.4041, df = 11.205, p-value = 0.1874
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -87.42307 397.42307
## sample estimates:
## mean of x mean of y 
##      1230      1075
t.test(wage_jp, wage_us, var.equal = TRUE)  # t test (assuming equal variance)
## 
##  Two Sample t-test
## 
## data:  wage_jp and wage_us
## t = 1.479, df = 16, p-value = 0.1586
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -67.16377 377.16377
## sample estimates:
## mean of x mean of y 
##      1230      1075

7.3.1.2 Titanic data

データ読み込み．

titanic <- read.csv("https://raw.githubusercontent.com/kurodaecon/bs/main/data/titanic3_csv.csv")

性別によって年齢の平均値が異なるかどうかの検定．

t.test(x = titanic$age[titanic$sex == "female"], y = titanic$age[titanic$sex == "male"])
## 
##  Welch Two Sample t-test
## 
## data:  titanic$age[titanic$sex == "female"] and titanic$age[titanic$sex == "male"]
## t = -2.0497, df = 798.36, p-value = 0.04072
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.71593356 -0.08034711
## sample estimates:
## mean of x mean of y 
##  28.68709  30.58523

t.test(変数 ~ グループ, data = dataset) で，2つのグループ間の検定が可能．

t.test(age ~ sex, data = titanic)  # same as above
## 
##  Welch Two Sample t-test
## 
## data:  age by sex
## t = -2.0497, df = 798.36, p-value = 0.04072
## alternative hypothesis: true difference in means between group female and group male is not equal to 0
## 95 percent confidence interval:
##  -3.71593356 -0.08034711
## sample estimates:
## mean in group female   mean in group male 
##             28.68709             30.58523

\(p < 0.05\) なので，「年齢の平均値は男性と女性で等しい」という帰無仮説を有意水準5%で棄却できる． 男性の標本平均の方が高いという事実と合わせると，「男性の方が有意に平均年齢が高い」と言える．

サブサンプルに対して検定する場合やデータフレームを加工して作成した変数を用いて検定する場合はパイプ演算子を使うと便利（tidyverse が必要）．

  • t.test の直前のパイプ演算子より前の部分（データフレーム）が t.test 関数の第一引数にならないため，t.test 関数内で data = . を明示的に書く必要がある（. がパイプ演算子の前の部分を受け取る）
library(tidyverse)
titanic %>% 
  filter(pclass == 1) %>%
  t.test(age ~ sex, data = .)
## 
##  Welch Two Sample t-test
## 
## data:  age by sex
## t = -2.3283, df = 278.85, p-value = 0.02061
## alternative hypothesis: true difference in means between group female and group male is not equal to 0
## 95 percent confidence interval:
##  -7.3664537 -0.6169013
## sample estimates:
## mean in group female   mean in group male 
##             37.03759             41.02927

7.3.2 Paired (dependent) \(t\) test

対応のある t 検定は t.test の引数で paired = TRUE 引数を追加して計算する．

wage_w <- c(1000, 1200, 1300, 1200, 1150, 1000, 1450, 1500, 1150, 1350)  # wife
wage_h <- c(900, 1300, 1200, 800, 1600, 850, 1000, 950, 1200, 1400)  # husband
# length(wage_w); length(wage_h)  # = 10
t.test(wage_w, wage_h, paired = TRUE)
## 
##  Paired t-test
## 
## data:  wage_w and wage_h
## t = 1.1638, df = 9, p-value = 0.2744
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -103.8108  323.8108
## sample estimates:
## mean difference 
##             110
t.test(wage_w - wage_h)  # same as an ordinary "one-sample t test"
## 
##  One Sample t-test
## 
## data:  wage_w - wage_h
## t = 1.1638, df = 9, p-value = 0.2744
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -103.8108  323.8108
## sample estimates:
## mean of x 
##       110

上記の通り「要素ごとの差」という一標本の検定になっている．

7.4 Multiple testing problem: Simulation / 多重検定の問題に関するシミュレーション

複数の標本（または複数の標本組み合わせ）に関する帰無仮説を独立して検定すると，そのうち少なくとも一つで誤って帰無仮説を棄却する確率が事前に設定した有意水準（5%など）を上回ってしまう．

1,000 の標本のうち「1つ」について \(H_0: \mu = 0\) を有意水準 \(\alpha = 0.05\) で検定した場合，帰無仮説が誤って棄却されてしまう理論上の割合は有意水準である5%．

\[ 1 - \mbox{(1つの標本について棄却されない確率)} = 1 - (1 - 0.05) = 0.05 \]

シミュレーション上は以下の通り（概ね5%）．

p_value <- NULL
for (i in 1:1000) {
  p_value_i <- t.test(rnorm(100, mean = 0, sd = 1))$p.value
  p_value <- c(p_value, p_value_i)
}
table(p_value <= 0.05)
## 
## FALSE  TRUE 
##   959    41

1,000 の標本のうち「2つ」について \(H_0: \mu = 0\) を検定した場合，少なくともどちらか一方で帰無仮説が誤って棄却されてしまう割合は以下の通り（概ね10%）．

\[ 1 - \mbox{(2回とも棄却されない確率)} = 1 - (1 - 0.05)^2 = 0.0975 \]

  • sapply(X = ベクトル, FUN = 関数) ：ベクトルの要素を1つずつ関数で評価しその結果をベクトルで返す
  • sample(ベクトル, 個数) ：指定した個数分だけベクトルからランダムにサンプリング
  • any(論理値ベクトル) ：論理値ベクトルの要素の1つ以上が TRUE の場合に TRUE を返す
table(sapply(X = 1:1000, FUN = function (x) any(sample(p_value, 2) <= 0.05)))
## 
## FALSE  TRUE 
##   918    82

1,000 の標本のうち「5つ」について \(H_0: \mu = 0\) を検定した場合，少なくとも1つで帰無仮説が誤って棄却されてしまう割合は以下の通り（概ね23%）．

\[ 1 - (1 - 0.05)^5 \approx 0.226 \]

table(sapply(X = 1:1000, FUN = function (x) any(sample(p_value, 5) <= 0.05)))
## 
## FALSE  TRUE 
##   796   204

7.4.1 ANOVA

同一母集団から3つの標本を作成し分散分析（ANOVA）を適用する．

df <- data.frame(val = rnorm(100 * 3), group = rep(1:3, 100))
head(df)
##          val group
## 1 -1.0417865     1
## 2  0.6511979     2
## 3 -0.1467194     3
## 4  0.1541756     1
## 5 -1.6550874     2
## 6  0.4495309     3
summary(aov(val ~ group, data = df))
##              Df Sum Sq Mean Sq F value Pr(>F)
## group         1   1.02  1.0234    1.07  0.302
## Residuals   298 284.98  0.9563
summary(aov(val ~ group, data = df))[[1]]$`Pr(>F)`[1]  # extract p-value
## [1] 0.3017551

以下のように \(p\) 値の分布を確認できる．

p_value <- NULL
for (i in 1:1000) {
  df_i <- data.frame(val = rnorm(100 * 3), group = rep(1:3, 100))
  p_i <- summary(aov(val ~ group, data = df_i))[[1]]$`Pr(>F)`[1]
  p_value <- c(p_value, p_i)
}
hist(p_value, breaks = 20, col = c("lightblue", rep("grey", 19)))

table(p_value <= 0.05)
## 
## FALSE  TRUE 
##   952    48

.