最大固有値のサンプリング分布のシミュレーション

Question

Xの最大固有値（Wishart分布に従う）を見つけたいと思います。そして、これらの固有値の経験的分布を見るためにシミュレーションを使用します。しかし、私はこのようにコード化するとき

library(MASS) 

function(X){ 
    maxeigen.XtX <- NULL 
    num_samples <- 1000 
    for(i in 1:num_samples){ 
     X <- mvrnorm(n=10,mu=rep(0,3),Sigma = matrix(c(1,0.2,0.1,0.2,1,0.2,0.1,0.2,1),nrow=3)) 
     XtX <- t(X)%*%X 
     maxeigen.XtX[i] <- max(eigen(XtX)$values) 
     } 
     return(maxeigen.XtX) 
     summary <- summary(maxeigen.XtX) 
     histgram <- hist(maxeigen.XtX,breaks=100) 
    } 

それは私の何かを与えるものではありません。どこに問題があるのか​​わからない？

Answer 1

はAがあなたの目標共分散行列とする：ここで

A <- matrix(c(1,0.2,0.1,0.2,1,0.2,0.1,0.2,1), nrow = 3) 

#  [,1] [,2] [,3] 
#[1,] 1.0 0.2 0.1 
#[2,] 0.2 1.0 0.2 
#[3,] 0.1 0.2 1.0 

は最大固有値のNサンプルを取得するための仕事関数です。 Aの行列分解は、N回ではなく1回だけ行われるので、MASS::mvrnormを使用するよりもはるかに効率的です。

g <- function (N, n, A) { 
    ## get upper triangular Choleksy factor of covariance `A` 
    R <- chol.default(A) 
    ## a function to generate `n` samples from `N(0, A)` 
    ## and get largest eigen value 
    f <- function (n, R) { 
    Xstd <- matrix(rnorm(n * dim(R)[1L]), n) ## `n` standard normal samples 
    X <- Xstd %*% R ## transform to have covariance `A` 
    S <- crossprod(X) ## `X'X` 
    max(eigen(S, symmetric = TRUE)$values) ## symmetric eigen decomposition 
    } 
    ## replicate `N` times for `N` samples of largest eigen values 
    replicate(N, f(n, R)) 
    } 

## try `N = 1000`, `n = 10`, as in your original code 
set.seed(0); x <- g(1000, 10, A) 

注：私はgに要約とプロットを依頼していません。サンプルがある限り、いつでも実行できます。

d <- density.default(x) ## density estimation 
h <- hist.default(x, plot = FALSE) ## histogram 
graphics:::plot.histogram(h, freq = FALSE, ylim = c(0, max(h$density, d$y)), 
          main = "histogram of largest Eigen value") 
lines(d$x, d$y, col = 2)