Моделирование распределения выборки самого большого собственного значения

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 образцов наибольшего значения собственному. Это намного эффективнее, чем использование MASS::mvrnorm, поскольку матричная факторизация A выполняется только один раз, а не N раз.

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)