0

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

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

生成相似的随机数据：

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% 的间隔。

4

2 回答 2

1

All you need is an interval that contains 95% of the generated b values. You can the highest posterior density from bayesian estimation, it's just that. There are many packages that compute it, for instance, function emp.hpd from TeachingDemos. Add

require(TeachingDemos)

and change the last line (return(b)) from ciint to

emp.hpd(b)

(No need to use return().)

于 2013-09-21T18:34:58.010 回答
0

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

于 2013-09-21T19:08:25.057 回答