bootパッケージを軽く触ってみたので、コードを残してみます。
ポアソン分布
パラメーター推定
ヘッシアンから信頼区間をつくると、正のパラメーターなのにマイナスが入ってしまう例を考えます。
set.seed(751) n <- 30 x <- rpois(n, lambda = 0.123) llf <- function(lambda){ sum(dpois(x, lambda, log = TRUE)) } r_optim <- optim(0.1, llf, method = "L-BFGS-B", lower = 1e-6, control = list(fnscale = -1), hessian = TRUE) se <- sqrt(diag(solve(-r_optim$hessian))) (r_optim_confint <- r_optim$par + se*qt(c(0.025, 0.975), df = n - 1))
[1] -0.03479004 0.10147670
漸近論から正規分布で近似できることになっているわけですが、サンプルサイズが足りないと困った値が出てきます。
ブートストラップで信頼区間を計算しましょう。ここではブートストラップ内で平均値を使って推定していますが、上のように最尤法を用いても構いません。
library(boot) r_boot <- boot(x, \(data, i){ mean(data[i]) }, R = 1000) (r_boot_confint <- quantile(r_boot$t[ ,1], c(0.025, 0.975)))
2.5% 97.5% 0.0 0.1
実際に取りうる値のサンプルから信頼区間をつくるので、正の値が出て来ます。
二標本ノンパラメトリック検定
多項分布から標本をつくります。
y <- rmultinom(n = 10, size = 1, prob = c(0.5, 0.25, 0.125, 0,125))
ポアソン分布の標本と同じ回数のブートストラッピングをかけます。
r_boot_mn <- boot(y, \(data, i){ mean(data[i]) }, R = 1000)
二標本の平均値の差を求めて集計します。
sum(r_boot$t[, 1] > r_boot_mn$t[, 1])/nrow(r_boot$t)
[1] 0.197
19.7%なので、ポアソン分布からの標本の平均値が、多項分布からの標本の平均値よりも大きいとも小さいとも言えないです。
線形回帰
ブートストラッピング向きではないのですが、線形回帰でも利用してみましょう。
信頼区間
データセットとモデルを設定します。
library(mvtnorm) df01 <- with(new.env(), { n <- 1000 mu <- c(4, 5) S <- matrix(c(2, -1, -1, 3), 2, 2) X <- rmvnorm(n, mu, S) colnames(X) <- c("x1", "x2") beta <- c(1, -1) y <- X %*% beta + rnorm(n, sd = 2) data.frame(y, X) }) nos <- 123 int <- c(0.025, 0.975) model <- y ~ x1 + x2
サンプルをそのまま使った線形回帰を行います。
r_lm <- lm(model, df01)
ブートストラッピングで係数の信頼区間を求めます。
library(boot) set.seed(418) r_boot <- boot(df01, function(data, i){ coef(lm(model, data[i, ])) }, nos) confint_B <- apply(r_boot$t, 2, quantile, prob = int) colnames(confint_B) <- names(coef(r_lm))
分散共分散行列から計算する一般的な信頼区間と、ブートストラッピング信頼区間を比較してみましょう。
(cmp_confint <- cbind(t(confint(r_lm)), confint_B))
(Intercept) x1 x2 (Intercept) x1 x2 2.5 % -0.9062976 0.9152777 -1.034447 -0.7924888 0.9116249 -1.0226879 97.5 % 0.4117309 1.1031148 -0.880196 0.5010901 1.0894776 -0.8808674
左が一般的な信頼区間、右がブートストラッピングです。
ブートストラッピングのサンプルサイズは123と大きくはないですが、ぼちぼち近いところになっています。
予測区間
予測区間も同様に出ますが、つくりが複雑なので手が込みます。賢い人がFAQサイトで説明してくれていたコードを踏襲しましょう。
# 予測する点 df02 <- data.frame(y = 0, x1 = 3, x2 = -1) # the value of y is not to be used. # 線型回帰で推定 prdct_lm <- predict(r_lm, df02[1, ], interval = "prediction") y_fit <- c(prdct_lm[, "fit"]) adjusted_residuals <- function(r_lm){ X <- model.matrix(r_lm) h <- X %*% solve(t(X) %*% X) %*% t(X) # 予測値の分散は diag(x^T σ^2 (X^t X)^{-1} x^T) = σ^2 {leverage} となるので、観測値と予測値と残>差を{leverage}で割って、観測値の分散σ^2に従う誤差のサンプルとする r <- residuals(r_lm) / sqrt(1 - diag(h)) # diag(h) is the leverage. r - mean(r) } adj_reg_lm <- adjusted_residuals(r_lm) set.seed(419) r_boot <- boot(adj_reg_lm, function(data, i){ df01$y <- y_fit + data[i] r_lm_boot <- lm(model, df01) adj_reg_boot <- adjusted_residuals(r_lm_boot) return( y_fit + # 非ブートストラッピング予測値 c(model.matrix(model, df02) %*% (coef(r_lm) - coef(r_lm_boot))) + # ブートストラッピングと 非〃の係数の推定量の差からの予測値の誤差 sample(adj_reg_boot, size = 1, replace = TRUE) # 観測値の分散に従う誤差 ) }, nos) prdct_B <- quantile(c(r_boot$t), prob = int)
データセットではなく、推定された誤差をリサンプリングしています。
比較すると、分散共分散行列との乖離はそこそこ大きいです。
(cmp_predict <- matrix(c(prdct_lm[2:3], prdct_B), 2, 2, dimnames = list(int, c("lm", "boot"))))
lm boot 0.025 0.3066421 0.6087671 0.975 8.2840626 7.8381543
このデータセットでは標準誤差が過小推定になっていますが、過大のときもあります。nosを増やしてどれぐらい改善するのかは、試してみてください。
