R中的自定义引导置信区间

Question

我需要找到一种方法来获取我使用自定义函数获得的估计值的引导置信区间。现在，问题是我有一个大矩阵，我从中随机取出行，然后计算所需的数量。

这是（希望）可重现的示例

生成相似的随机数据：

mat1 <- matrix(rnorm(300, 80, 20), nrow = 100)

计算所需量的函数（其中 R 是相关矩阵）：

IIvar <- function(R) { 
d <- eigen(R)$values  
p <- length(d)  
sum((d-1)^2)/(p*(p-1))}

我尝试求解的函数（其中 omat 是由一些 mat1 行组成的较小矩阵，freq 是 omat 中的行数，numR 是复制数）：

ciint <- function(omat, mat1, freq, numR) {
II <- IIvar(cor(omat))
n <- dim(mat1)[1]
b <- numeric(numR)
for (i in 1:numR) { b[i] <- IIvar(cor(mat1[sample(c(1:n),freq),]))}
hist(b)
abline(v = II, lty = 5, lwd = 3)
return(b) }

结果向量 b 具有从 mat1 中随机选择的行（数量由 freq 确定）的矩阵获得的所有值，可以与来自 omat 的 IIvar（由总体成员资格选择的行的矩阵）进行比较。

在 mat1 中，我有 5 个群体（按行分组），我需要分别计算所有这些群体的 IIvar 并为获得的值生成置信区间。

当我像这样运行我的 ciint 函数时

ciint(omat, mat1, 61, 1000)

我得到了值的分布，以及“真实”IIvar 值的位置，但我不知道如何从这一点生成 95% 的区间。

Answer 1

您只需要一个包含 95% 生成的b 值的区间。你可以从贝叶斯估计中得到最高的后验密度，就是这样。有许多计算它的包，例如，来自TeachingDemos 的函数emp.hpd。添加

require(TeachingDemos)

并将最后一行 (return(b)) 从 ciint 更改为

emp.hpd(b)

（无需使用return()。）

Answer 2

我不确定你想用你的函数来完成什么，但如果你想做 boostrapping，那么看看boot 包中的boot 函数。您可以将自定义函数传递给boot，它将获取引导样本，将它们传递给自定义函数，然后整理结果。然后，它还具有多个结果的置信区间选项。