Le 16/02/2012 16:20, [email protected] a écrit :
I don't see any way to fix multivariate_normal for this case, except
for dropping svd or for random perturbing a covariance matrix with
multiplicity of singular values.
Hi,
I just made a quick search in what R guys are doing. It happens there
are several codes (http://cran.r-project.org/web/views/Multivariate.html
). For instance, mvtnorm
(http://cran.r-project.org/web/packages/mvtnorm/index.html). I've
attached the related function from the source code of this package.
Interestingly enough, it seems they provide 3 different methods (svd,
eigen values, and Cholesky).
I don't have the time now to dive in the assessments of pros and cons of
those three. Maybe one works for our problem, but I didn't check yet.
Pierre
# $Id: mvnorm.R 222 2011-01-31 14:02:02Z thothorn $
rmvnorm<-function (n, mean = rep(0, nrow(sigma)), sigma = diag(length(mean)),
method=c("eigen", "svd", "chol"))
{
if (!isSymmetric(sigma, tol = sqrt(.Machine$double.eps),
check.attributes = FALSE)) {
stop("sigma must be a symmetric matrix")
}
if (length(mean) != nrow(sigma)) {
stop("mean and sigma have non-conforming size")
}
sigma1 <- sigma
dimnames(sigma1) <- NULL
if(!isTRUE(all.equal(sigma1, t(sigma1)))){
warning("sigma is numerically not symmetric")
}
method <- match.arg(method)
if(method == "eigen"){
ev <- eigen(sigma, symmetric = TRUE)
if (!all(ev$values >= -sqrt(.Machine$double.eps) * abs(ev$values[1]))){
warning("sigma is numerically not positive definite")
}
retval <- ev$vectors %*% diag(sqrt(ev$values),
length(ev$values)) %*% t(ev$vectors)
}
else if(method == "svd"){
sigsvd <- svd(sigma)
if (!all(sigsvd$d >= -sqrt(.Machine$double.eps) * abs(sigsvd$d[1]))){
warning("sigma is numerically not positive definite")
}
retval <- t(sigsvd$v %*% (t(sigsvd$u) * sqrt(sigsvd$d)))
}
else if(method == "chol"){
retval <- chol(sigma, pivot = TRUE)
o <- order(attr(retval, "pivot"))
retval <- retval[,o]
}
retval <- matrix(rnorm(n * ncol(sigma)), nrow = n) %*% retval
retval <- sweep(retval, 2, mean, "+")
colnames(retval) <- names(mean)
retval
}
dmvnorm <- function (x, mean, sigma, log=FALSE)
{
if (is.vector(x)) {
x <- matrix(x, ncol = length(x))
}
if (missing(mean)) {
mean <- rep(0, length = ncol(x))
}
if (missing(sigma)) {
sigma <- diag(ncol(x))
}
if (NCOL(x) != NCOL(sigma)) {
stop("x and sigma have non-conforming size")
}
if (!isSymmetric(sigma, tol = sqrt(.Machine$double.eps),
check.attributes = FALSE)) {
stop("sigma must be a symmetric matrix")
}
if (length(mean) != NROW(sigma)) {
stop("mean and sigma have non-conforming size")
}
distval <- mahalanobis(x, center = mean, cov = sigma)
logdet <- sum(log(eigen(sigma, symmetric=TRUE,
only.values=TRUE)$values))
logretval <- -(ncol(x)*log(2*pi) + logdet + distval)/2
if(log) return(logretval)
exp(logretval)
}
# file MASS/R/mvrnorm.R
# copyright (C) 1994-2004 W. N. Venables and B. D. Ripley
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 or 3 of the License
# (at your option).
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# A copy of the GNU General Public License is available at
# http://www.r-project.org/Licenses/
#
mvrnorm <- function(n = 1, mu, Sigma, tol=1e-6, empirical = FALSE)
{
p <- length(mu)
if(!all(dim(Sigma) == c(p,p))) stop("incompatible arguments")
eS <- eigen(Sigma, symmetric = TRUE, EISPACK = 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) # remove means
X <- X %*% svd(X, nu = 0)$v # rotate to PCs
X <- scale(X, FALSE, TRUE) # rescale PCs to unit variance
}
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)
}
_______________________________________________
NumPy-Discussion mailing list
[email protected]
http://mail.scipy.org/mailman/listinfo/numpy-discussion
