R中置信区间的模拟

Question

我有一个 R 函数，它为 t 分布的 ncp（非中心性参数）提供 95% 置信区间。

通过 R 中的模拟，是否有可能表明，从长远来看，来自该 R 函数的 CI 捕获给定的 TRUE ncp（此处的“2”与输入 t) 95% 的时间？

（我很欣赏关于如何做到这一点的任何想法）

CI.ncp <- function(t, N){

  f <- function (ncp, alpha, q, df) {  
abs(suppressWarnings(pt(q = t, df = N - 1, ncp, lower.tail = FALSE)) - alpha) }

sapply(c(0.025, 0.975),
function(x) optim(1, f, alpha = x, q = t, df = N - 1, control = list(reltol = (.Machine$double.eps)))[[1]]) }

#Example of Use:
CI.ncp(t = 2, N = 20) # gives: -0.08293755  4.03548862 

#(in the long-run 95% of the time, "2" is contained within these
# two numbers, how to show this in R?)

这是我尝试过但没有成功的方法：

fun <- function(t = 2, N = 20){

  ncp = rt(1, N - 1, t)
  CI.ncp(t = 2, N = 20)
  mean(ncp <= 2 & 2 <= ncp )
   }

 R <- 1000
 sim <- t(replicate(R, fun()))
 coverage <- mean(sim[,1] <= 2 & 2 <= sim[,2])

Answer 1

问题是我们需要在CI.ncp 中提供从fun 获得的随机ncp：

 fun <- function(t = 2, N = 20){ ;

   ncp = rt(1, N - 1, t);
   CI.ncp(t = ncp, N = 20);
   }

  R <- 1e4 ;
  sim <- t(replicate(R, fun()));
  coverage <- mean(sim[,1] <= 2 & 2 <= sim[,2])

Answer 2

我会使用包MBESS。

#install.packages("MBESS")
library(MBESS)

fun <- function(t = 2, N = 20, alpha = 0.95){
   x = rt(1, N - 1, t)
   conf.limits.nct(x, df = N, conf.level = alpha)[c(1, 3)]
}

set.seed(5221)
R <- 1000
sim <- t(replicate(R, fun()))
head(sim)
coverage <- mean(sim[,1] <= 2 & 2 <= sim[,2])
coverage
[1] 0.941