5

私の意図は、2 つの共分散行列間の全体的な類似性を見つけることを目的としたいくつかの関数を作成することでした。これには、それらをランダムなベクトルで乗算して応答ベクトルを相関させるか、行列の 1 つをブートストラップして比較に役立つ相関係数分布を取得することによって行います。しかし、どちらの場合も、誤った結果が得られました。観測されたマトリックス間の相関は 0.93 まで高かったが、分布は最大 0.2 までしか範囲がなかった。これは関数のコードです:

resamplerSimAlt <- function(mat1, mat2, numR, graph = FALSE)
{
  statSim <- numeric(numR)
  mat1vcv <- cov(mat1)
  mat2vcvT <- cov(mat2)
  ltM1 <- mat1vcv[col(mat1vcv) <= row(mat1vcv)]
  ltM2T <- mat2vcvT[col(mat2vcvT) <= row(mat2vcvT)]
  statObs <- cor(ltM1, ltM2T)                           
  indice <- c(1:length(mat2))
  resamplesIndices <- lapply(1:numR, function(i) sample(indice, replace = F))
  for (i in 1:numR)
  {
    ss <- mat2[sample(resamplesIndices[[i]])]
    ss <- matrix(ss, nrow = dim(mat2)[[1]], ncol = dim(mat2)[[2]])
    mat2ss <- cov(ss)
    ltM2ss <- mat2ss[col(mat2ss) <= row(mat2ss)]
    statSim[i] <- cor(ltM1, ltM2ss)
  }
  if (graph == TRUE)
  {
    plot(1, main = "resampled data density distribution", xlim = c(0, statObs+0.1), ylim = c(0,14))
    points(density(statSim), type="l", lwd=2)
    abline(v = statObs)
    text(10, 10, "observed corelation = ")
  }
  list( obs = statObs , sumFit = sum(statSim > statObs)/numR)
}

実際、2 つの元のマトリックス間の相関係数が高く、1 番目の元のマトリックスと 2 番目の再サンプリングされたマトリックスの間の相関係数が、10000 回のブートストラップ反復後に最大 0.2 であるとは信じがたいです。

コードの有効性に関するコメントはありますか?

4

1 に答える 1

2

申し訳ありませんが、私は 2 つの共分散行列間の効率的な相関をチェックするという目標に追いつくのに十分な教育を受けていませんが、コード自体を理解しようとしました。

私が正しければ、同じ行列 ( ) から 10.000 の異なる行列を作成していることになります。セル全体を並べ替え、の共分散行列とリサンプリングされた配列の共mat2分散行列の間の相関を再計算することです。mat1それらはstatSim変数に格納されます。

元の相関効率が高かった ( statObs) とおっしゃいましたが、最大値statSimが低かったのは奇妙です。問題は結果リストにあると思います：

list( obs = statObs , sumFit = sum(statSim > statObs)/numR)

元の相関係数 ( ) を返しますobsが、 で書かれた最大値ではありませんsumFit。そこでは、たとえば使用できます。max(statSim). リサンプリングによって相関が効率的に改善されたかどうかを確認するために戻るポイントはわかりsumFitますが、コードに基づいて、理論に問題はありません。

maxシミュレートされた相関係数の更新された関数:

resamplerSimAlt <- function(mat1, mat2, numR, graph = FALSE)
{
  statSim <- numeric(numR)
  mat1vcv <- cov(mat1)
  mat2vcvT <- cov(mat2)
  ltM1 <- mat1vcv[col(mat1vcv) <= row(mat1vcv)]
  ltM2T <- mat2vcvT[col(mat2vcvT) <= row(mat2vcvT)]
  statObs <- cor(ltM1, ltM2T)                           
  indice <- c(1:length(mat2))
  resamplesIndices <- lapply(1:numR, function(i) sample(indice, replace = F))
  for (i in 1:numR)
  {
    ss <- mat2[sample(resamplesIndices[[i]])]
    ss <- matrix(ss, nrow = dim(mat2)[[1]], ncol = dim(mat2)[[2]])
    mat2ss <- cov(ss)
    ltM2ss <- mat2ss[col(mat2ss) <= row(mat2ss)]
    statSim[i] <- cor(ltM1, ltM2ss)
  }
  if (graph == TRUE)
  {
    plot(1, main = "resampled data density distribution", xlim = c(0, statObs+0.1), ylim = c(0,14))
    points(density(statSim), type="l", lwd=2)
    abline(v = statObs)
    text(10, 10, "observed corelation = ")
  }
  list( obs = statObs , sumFit = sum(statSim > statObs)/numR, max=max(statSim))
}

私が実行したもの：

> mat1 <- matrix(runif(25),5,5)
> mat2 <- mat1+0.2
> resamplerSimAlt(mat1, mat2, 10000)
$obs
[1] 1

$sumFit
[1] 0

$max
[1] 0.94463

そしてランダムでmat2：

> mat2 <- matrix(runif(25),5,5)
> resamplerSimAlt(mat1, mat2, 10000)
$obs
[1] 0.31144

$sumFit
[1] 0.9124

$max
[1] 0.9231

私の答えは本当の答えではないかもしれません。その場合は、問題の詳細をお知らせください:)

于 2011-10-23T15:19:08.333 に答える