【R】スマートに検定を使う前に数値実験で考えてみる - ケィオスの時系列解析メモランダム

今回は、理論的に導かれる検定手法を適用する前に、数値実験（モンテカルロ・シミュレーション）を通して、検定の考え方を考察します。

母比率の推定

　母比率（population proportion）とは、母集団の中で、ある条件を満たす要素の割合のことです。

　例えば、

日本全国の成人のうち、毎日コーヒーを飲む人の割合
アンケート回答者の中で「満足」と答えた人の割合

などが母比率です。

記号で表すと、母比率は $\displaystyle{ p }$ で表されます。

　ここでは、標本（sample）を取り、そのデータから計算された比率で母比率を推定します。この標本の比率を 標本比率（sample proportion）といいます。

　標本比率は以下の式で求めます：

$\displaystyle{ \hat{p} = \frac{x}{n} }$

$\displaystyle{ x }$ ：標本の中で条件を満たす数
$\displaystyle{ n }$ ：標本の大きさ（標本の人数など）

　標本比率 $\displaystyle{ \hat{p} }$ は、標本の違いで当然ばらつきます。

　どの程度ばらつくかを調べるために、以下のスクリプトを試してください。


# ───────────────────────────────────────────────────────────────
# 母比率の推定：標本から推定した比率の分布を描くスクリプト
# ───────────────────────────────────────────────────────────────

# 1. パラメータ設定 ------------------------------------------------
p      <- 0.3     # 母比率（真の成功確率）
N      <- 50      # 標本サイズ
n.rep  <- 1000   # 繰り返し試行回数

# 2. 推定比率格納用ベクトルの初期化 --------------------------------
phat.rep <- numeric(n.rep)

# 3. シミュレーションループ ------------------------------------------
for (i in seq_len(n.rep)) {
  # Bernoulli試行を N 回行い、成功の割合を計算
  x            <- rbinom(N, size = 1, prob = p)
  phat.rep[i]  <- mean(x)
}

# 4. 描画範囲の設定 ------------------------------------------------
# 理論的な正規近似のための平均と標準偏差
mean_phat <- p
sd_phat   <- sqrt(p * (1 - p) / N)

# 5. ヒストグラムの描画 --------------------------------------------
par(mfrow = c(1, 1), mar = c(4, 4, 2, 1))
hist(phat.rep,
     breaks = seq(0, 1, by = 0.05),  # ビン幅0.01
     freq   = FALSE,                 # 密度表示
     col    = "lightblue",
     xlab   = "Sample proportion (\\hat{p})",
     xlim   = c(0, 1),
     main   = sprintf("Sampling distribution of \\hat{p} (p=%.2f, N=%d)", p, N))

# 6. 実証的密度曲線を重ねる ----------------------------------------
lines(density(phat.rep), lwd = 2, col = "blue", lty = 1)

# 7. 理論的正規近似を重ねる ----------------------------------------
x_vals <- seq(0, 1, length.out = 1000)
curve(dnorm(x, mean = mean_phat, sd = sd_phat),
      from = 0, to = 1,
      add  = TRUE, lwd = 2, col = "red", lty = 2)

# 8. 凡例の追加 ----------------------------------------------------
legend("topright",
       legend = c("Empirical density", "Normal approx."),
       col    = c("blue",      "red"),
       lwd    = 2,
       lty    = c(1, 2),
       bty    = "n")

　この数値実験を通じて、もし真の $p$ (母比率)を知っていたとして、標本数が $N$ のもとで得られる標本比率のばらつきの様子を確認することができます。下の図の上段は、その理想的なケースを示しています。

　一方、実際の推定においては、真の $p$ 未知です。そのため、推定された標本比率 $\hat{p}$ から、真の $p$ がどれだけ離れているかを考える必要があります。図の下段のように、標本比率のばらつきの大きさがわかれば、母比率の推定誤差の大きさについて推論することが可能になります。

z統計量とt統計量の分布

　平均値の推定では、z統計量とt統計量という指標が登場します。これらは、標本平均が母平均からどの程度離れているかを標準化して評価するための統計量です。

z統計量（z-score）

z統計量は、母分散（あるいは母標準偏差）が既知である場合に使用されます。標本平均 $\displaystyle{ \bar{x} }$ の分布が正規分布に従うと仮定し、以下の式で計算されます：

$\displaystyle{ z=\frac{\bar{x}-\mu}{\sigma / \sqrt{n}} }$

ここで、 $\displaystyle{ \bar{x} }$ は標本平均、 $\mu$ は母平均、 $\sigma$ は母標準偏差 (真の標準偏差を既知としている)、 $n$ は標本サイズです。

　この z 統計量は標準正規分布（平均0、分散1）に従うとされ、母平均に関する検定や信頼区間の推定に使われます。

t統計量（t-score）

　一方、現実の多くの場面では母標準偏差 $\sigma$ が不明です。事前に私たちがすることができるのは標本サイズ $n$ だけです。

　その場合は標本から計算された不偏分散の平方根

$\displaystyle{ s=\sqrt{\frac{1}{n-1} \sum_{i=1}^n\left(x_i-\bar{x}\right)^2} }$

を用いて、t統計量を使用します：

$\displaystyle{ t=\frac{\bar{x}-\mu}{s / \sqrt{n}} }$

　この t統計量は、自由度 $n−1$ の t分布に従います。t分布は、標本サイズが小さいときには裾が厚い（極端な値が出やすい）という特徴がありますが、標本サイズが大きくなるにつれて標準正規分布に近づきます。

　実際には、真の標準偏差 $\sigma$ がどの程度かは、ほとんどの場合わかりません。そのため、標本から得られた標準偏差の推定値には誤差が含まれ、その不確かさを統計解析において考慮する必要があります。

　特に標本数が少ないときには、推定された標準偏差 $s$ が真の値 $\sigma$ よりも過小評価されることがあり、その結果、標本平均の推定結果の誤差も大きくなる可能性があります。したがって、信頼区間の幅や検定の判断において、t分布を用いることでこの不確かさを適切に補正することが求められます。


# ------------------------------------------------------------
# t統計量の分布：母分散既知 vs 母分散未知（breaks 揃え・コメント付き）
# ------------------------------------------------------------

# 1. パラメータ設定 ------------------------------------------
n_sample <- 6           # 各試行における標本サイズ
n_trials <- 10000       # 試行回数（シミュレーションの繰り返し数）
mu <- 0                 # 母平均（真の平均値）
sigma <- 2              # 母標準偏差（既知と仮定）

# 2. 結果保存用ベクトル初期化 ------------------------------
z_stats <- numeric(n_trials)  # 母分散既知の場合の統計量（z）
t_stats <- numeric(n_trials)  # 母分散未知の場合の統計量（t）

# 3. シミュレーションによる統計量の生成 ---------------------
for (i in 1:n_trials) {
  x <- rnorm(n_sample, mean = mu, sd = sigma)  # 標本を正規分布から抽出
  x_bar <- mean(x)                             # 標本平均
  s <- sd(x)                                   # 標本標準偏差

  # z統計量（母分散既知）
  z_stats[i] <- (x_bar - mu) / (sigma / sqrt(n_sample))

  # t統計量（母分散未知）
  t_stats[i] <- (x_bar - mu) / (s / sqrt(n_sample))
}

# 4. 共通ビン（breaks）設定 ---------------------------------
# ヒストグラムのビンを統一することで、両者の分布を比較しやすくする
xlim_range <- c(min(z_stats, t_stats), max(z_stats, t_stats))  # 範囲
n_breaks <- 200
breaks_shared <- seq(xlim_range[1], xlim_range[2], length.out = n_breaks + 1)

# 5. ヒストグラムと理論分布の描画 ----------------------------
par(mfrow = c(2, 1))  # 上下2つのグラフに分割表示

## （1）z統計量の分布（母分散既知）--------------------------
hist(z_stats, breaks = breaks_shared, col = rgb(1, 0, 0, 0.5),
     freq = FALSE, xlim = c(-4, 4),
     ylim = c(0, dnorm(0) * 1.2),
     main = paste0("z統計量の分布（母分散既知, n = ", n_sample, "）"),
     xlab = "統計量の値")

curve(dnorm(x), col = "red", lwd = 2, lty = 2, add = TRUE)  # 理論：標準正規分布

legend("topright",
       legend = c("シミュレーション", "標準正規分布"),
       fill = c(rgb(1, 0, 0, 0.5), NA), border = NA,
       lty = c(NA, 2), col = c(NA, "red"), lwd = c(NA, 2),
       bg = "white")

## （2）t統計量の分布（母分散未知）---------------------------
hist(t_stats, breaks = breaks_shared, col = rgb(0, 0, 1, 0.5),
     freq = FALSE, xlim = c(-4, 4),
     ylim = c(0, dt(0, df = n_sample - 1) * 1.2),
     main = paste0("t統計量の分布（母分散未知, n = ", n_sample, "）"),
     xlab = "統計量の値")

curve(dt(x, df = n_sample - 1), col = "black", lwd = 2, add = TRUE)  # 理論：t分布
curve(dnorm(x), col = "red", lwd = 2, lty = 2, add = TRUE)           # 理論：正規分布（比較用）

legend("topright",
       legend = c("シミュレーション", 
                  paste0("t分布（df=", n_sample - 1, "）"), 
                  "標準正規分布"),
       fill = c(rgb(0, 0, 1, 0.5), NA, NA), border = NA,
       lty = c(NA, 1, 2), col = c(NA, "black", "red"), lwd = c(NA, 2, 2),
       bg = "white")