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

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

这里是（希望）再现的示例

生成类似随机数据：

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))}

My功能（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随机选择的行（由数频率确定的）的矩阵可与IIvar从OMAT（与由人口成员选择的行的矩阵）进行比较得到的所有值。

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

当我像这样运行

ciint(omat, mat1, 61, 1000)

我得到的值的分布，而“真正的” IIvar值的位置我ciint功能，但我不知道如何从这个95周％的时间间隔点。

来源

2013-09-21 Ian Stuart

所有你需要的是一个包含95％产生的b值的区间。贝叶斯估计可以得到最高的后验密度，就是这样。有很多计算它的软件包，例如，从TeachingDemos起的功能emp.hpd。添加

require(TeachingDemos)

，并从ciint最后一行（return(b)）更改为

emp.hpd(b)

（无需使用return()）。

来源

2013-09-21 18:34:58

伟大的建议，这正是我需要的。与此同时，我发现这个[网站]（http://codealamode.blogspot.com/2013/08/bootstrap-confidence-interval-methods.html），其中列出了替代bootrstrap配置项以及R代码。你是否知道另一个与你的建议类似的软件包？ –

我不知道你想与完成什么你的功能，但如果你想做增压包装，那么看看boot包中的boot功能。您可以将自定义函数传递到boot，它将引导自举样本，将它们传递给自定义函数，然后整理结果。从结果中，它也有多个置信区间的选项。

来源

2013-09-21 19:08:25

引导包很棒，但我觉得它的语法有些令人费解。我无法让它从一个更大的矩阵中选择随机子集并计算出我需要的东西。 –

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

回答

相关问题