0

カスタム関数で取得した推定値のブートストラップ信頼区間を取得する方法を見つける必要があります。さて、問題は、ランダムに行を取り出して必要な量を計算する 1 つの大きな行列があることです。

これが（うまくいけば）再現可能な例です

同様のランダム データを生成します。

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

必要なのは、生成されたb値の 95% を含む間隔だけです。ベイジアン推定から最高の事後密度を得ることができます。それだけです。それを計算する多くのパッケージがあります。たとえば、関数emp.hpdfromTeachingDemosです。追加

require(TeachingDemos)

最後の行 ( return(b)) を からciintに変更します。

emp.hpd(b)

(使用する必要はありませんreturn()。)

于 2013-09-21T18:34:58.010 に答える