10年以上前に、「信頼区間の意味と、Rのpredict()関数の使い方の注意点」というエントリを書いてましたが、余計な記述が多かったのであらためて論点を整理しておきます。
統計学の教科書をまじめに読んでいる人にとっては今さらな話なのですが、「95%信頼区間」というのは、「この範囲に真の値（母数）がある確率が95%」ということを表しているのではありません。ところが、学生のような初学者だけでなく、理系の大学教員でも「この範囲に母数が含まれる確率が95%」というふうに間違った説明をしてしまっていることがあるので、意外と厄介な概念です。
たとえば単回帰分析を行って、約2.5という回帰係数が得られ、その95%信頼区間が約1.8〜3.2だったとします。
ここで、「真のパラメータが1.8〜3.2である確率が95%」と言ってしまうと間違いだということになります。頻度主義の統計学は、「真の値が所与の場合に、データがどう変動しそうか」を考えるものなので、「真の値がこの範囲で動く」という考え方は基本的には持ち出しません。ベイズ統計学の場合は、「データが所与の場合に、真の値はどのあたりにあると考えるのが確からしいか」と考えるので、真の値が動きますが、ベイズ統計学の立場を取っているわけではない文脈でというか、おそらくベイズ統計学を意識的に勉強したわけでもないのに、そういうイメージで語ってしまっている人がたまによく居ます。
頻度主義における「信頼区間」を正しくイメージするなら、
- 真のパラメータ（母数）を設定してやると、それに従っサンプルデータを生成する、乱数発生器を考える。
- サンプルを生成するたびに、一定の手続きに従って、そのサンプル（データセット）から母数を推定し、信頼区間も計算する。
- この作業を何回も繰り返すと、信頼区間そのものは、毎回異なるものが得られる。「サンプリングと信頼区間計算」を100回繰り返すと、100個の信頼区間が得られる。
- この100個の信頼区間のうち、95個ぐらいは、最初に設定した真のパラメータを含んでいる。
ということになります。一番のポイントは、手元に得られたデータから計算された1回1回の信頼区間の値そのものには、あまり意味はないということです。「今回得られた区間に95%収まるように何かが動く」と考えるから、間違えるわけです。そうではなく、サンプリングを繰り返すたびに「信頼区間自体が動く」んですよね。
別の言い方をすると、「95%」というのはなんとなく強そうな数値ですが、「今回得られた信頼区間」が95%の強さを持っているわけではないんです。そうではなく、「信頼区間を計算する手続き」（ルール）が、何度も実行されたときに、95%の確率で真の値を捕捉することができるという強さを持っているということです。95%というのは、いわば「信頼区間算出ルールの勝率」です。この、「算出ルールの勝率」という概念を持っておくと、うまく頭が整理されるのではないかと思います。
まぁ、「真のパラメータが1.8〜3.2である確率」を意味するベイズ信用区間も、標準的な設定（共役事前分布を使うとか）で計算すると結果的には一致するし、そっちのほうが直観的に理解しやすいのですが、頻度主義の統計手法を使っている限りは、明確に異なる概念であるということは理解しておく必要があります。
じつは上述の信頼区間は、が真のパラメータだとして、
という帰無仮説から
という帰無仮説までは、今回のデータによって棄却されないという意味にもなります。また、真のパラメータが1.8や3.2のとき、そのモデルからサンプリングを繰り返してパラメータを推定する操作を繰り返すとした場合に、真のパラメータからサンプリングされたデータから計算されるパラメータ推定値のばらつきの両側95%の範囲内に、今回の手元のデータから推定されたパラメータが入っているということでもあります（通常の線形回帰の場合）。
以前のエントリではその説明を先にしていたのですが、この言い方だと、真の値が1.8から3.2までの間で動くかのような印象を与えるので、やめといたほうがいいですね。この説明においても、あくまで「真のパラメータが1.8に固定されている場合」から「真のパラメータが3.2に固定されている場合」までを考えているということであって、真のパラメータを変動するものとして考えているわけではないのですが、ややこしいです（笑）
上の図を描いたコードは下記のとおり。回帰係数と信頼区間が特定の値になるように無理やり選んでいます。
library(ggplot2) set.seed(1) # 設定 target_beta <- 2.5 target_lwr <- 1.8 target_upr <- 3.2 target_half <- (target_upr - target_lwr) / 2 # 0.7 tol_beta <- 0.002 # 傾きの許容誤差 tol_half <- 0.002 # 半幅の許容誤差 max_tries <- 20000 # 試行回数 # x と ノイズ規模 n <- 50 x <- seq(1, 25, length.out = n) beta0 <- 0 beta1_true <- 2.5 half_width <- target_half # 0.7 Sxx <- sum((x - mean(x))^2) tcrit <- qt(0.975, df = n - 2) # sigma を逆算 sigma <- half_width * sqrt(Sxx) / tcrit # 探索ループ（中心と半幅で評価） best <- list(score = Inf) for (i in 1:max_tries) { y <- beta0 + beta1_true * x + rnorm(n, 0, sigma) df <- data.frame(x = x, y = y) fit <- lm(y ~ x, data = df) beta_hat <- unname(coef(fit)["x"]) ci <- confint(fit)["x", ] lwr <- unname(ci[1]) upr <- unname(ci[2]) half <- (upr - lwr) / 2 # 距離（中心と半幅だけで評価） score <- abs(beta_hat - target_beta) + abs(half - target_half) # ベスト更新 if (score < best$score) { best <- list( score = score, df = df, fit = fit, beta_hat = beta_hat, lwr = lwr, upr = upr, half = half, iter = i ) } # 十分近ければ終了 if (abs(beta_hat - target_beta) <= tol_beta && abs(half - target_half) <= tol_half) { message("Found at iter = ", i) break } } # 結果表示 best$iter best$beta_hat c(best$lwr, best$upr) best$half # 予測値＋信頼区間（塗りつぶし） newx <- data.frame(x = seq(min(best$df$x), max(best$df$x), length.out = 200)) pred <- predict(best$fit, newdata = newx, interval = "confidence", level = 0.95) pred_df <- cbind(newx, as.data.frame(pred)) # fit/lwr/upr # プロット ggplot() + geom_point(data = best$df, aes(x = x, y = y)) + geom_ribbon( data = pred_df, aes(x = x, ymin = lwr, ymax = upr), alpha = 0.25, inherit.aes = FALSE ) + geom_line( data = pred_df, aes(x = x, y = fit), linewidth = 1, inherit.aes = FALSE ) + labs( title = paste0( "Slope = ", round(best$beta_hat, 4), " / 95% CI = [", round(best$lwr, 4), ", ", round(best$upr, 4), "]" ), x = "x", y = "y" )
（参考：このブログのおすすめ記事一覧はコチラ）
