首页 文章

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

提问于
浏览
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%的间隔 .

2 回答

  • 1

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

    require(TeachingDemos)
    

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

    emp.hpd(b)
    

    (无需使用 return() . )

  • 0

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

相关问题