最大特征值采样分布的模拟

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)

最大特征值采样分布的模拟

Simulation of the sampling distribution of the largest eigen value

simulation

r

function

eigenvalue