首页 > 解决方案 > MASS 包中 MVRNorm 背后的代数

问题描述

这是代码:

function (n = 1, mu, Sigma, tol = 1e-06, empirical = FALSE, EISPACK = FALSE) 
{
    p <- length(mu)
    if (!all(dim(Sigma) == c(p, p))) 
        stop("incompatible arguments")
    if (EISPACK) 
        stop("'EISPACK' is no longer supported by R", domain = NA)
    eS <- eigen(Sigma, symmetric = TRUE)
    ev <- eS$values
    if (!all(ev >= -tol * abs(ev[1L]))) 
        stop("'Sigma' is not positive definite")
    X <- matrix(rnorm(p * n), n)
    if (empirical) {
        X <- scale(X, TRUE, FALSE)
        X <- X %*% svd(X, nu = 0)$v
        X <- scale(X, FALSE, TRUE)
    }
    X <- drop(mu) + eS$vectors %*% diag(sqrt(pmax(ev, 0)), p) %*% 
        t(X)
    nm <- names(mu)
    if (is.null(nm) && !is.null(dn <- dimnames(Sigma))) 
        nm <- dn[[1L]]
    dimnames(X) <- list(nm, NULL)
    if (n == 1) 
        drop(X)
    else t(X)
}

我很好奇的问题是:

 x <- eS$vectors %*% diag(sqrt(ev)) %*% t(x) # ignoring drop(mu)
 ...
 t(x)

为什么会这样

X^T = UVZ^T,其中 Z 是标准化的 MVN?

我原以为这将是 X = UVZ,其中 X ~ MVN(0, UV(I)(UV)^T) = MVN(0, Sigma)?

回应 Siong Thye Goh 的回答:

我可以看到代数,并且仅通过考虑维度就可以做到这一点,但是考虑到多元正态的属性,转置一切的整个行为似乎很奇怪。即 X = UVZ

我做了一些审查,发现这实际上是一个Matrix Normal,并且那里的仿射变换以类似的方式工作。也就是说,X = Z (UV)^T。

我不确定我在理解这一点时是否缺少一些愚蠢的东西,或者我是否完全错过了关于为什么所有内容都被转置的图片,例如MVN 的 Wikipedias 仿射变换

标签: r

解决方案


U 是 Sigma 的特征向量。即 Sigma = UV^2 U^T,其中 V 是对角矩阵。

让我们计算协方差矩阵 E[X^TX],看看它是否等于 Sigma,其中 X=UVZ^T 和 Z^T 满足 E[Z^TZ]=I,即单位矩阵。

我们有

E[X^TX]=E[UVZ^TZVU^T]=UVE[Z^TZ]VU^T=UV^2U^T=Sigma


推荐阅读