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中的自定義引導置信區間

回答

相關問題