R devel - Mar 2012 - Conditional means and variances of a multivariate normal distribution

Hi,



Let X = (x_1, x_2, ... ,  x_p) be multivariate normal with mean, mu = (mu_1, ...
, mu_p) and covariance = Sigma.  I was looking for an R function to compute
conditional mean and conditional variance of a given subset of X given another
subset of X.  While this is trivially easy to do, there is nothing in
"base" for doing this, at least nothing that I am aware of.  I am also
not aware of anything in the contributed packages (although my search was not
comprehensive).  I feel that this would be a useful addition, if it is not
already there.  I have written this following function, which I am sure can be
improved a lot (including better argument names!). I would like to hear your
thought on this.



condNormal <- function(x.given, mu, sigma, req.ind, given.ind){

# Returns conditional mean and variance of x[req.ind]

# Given x[given.ind] = x.given

# where X is multivariate Normal with

# mean = mu and covariance = sigma

#

B <- sigma[req.ind, req.ind]

C <- sigma[req.ind, given.ind]

D <- sigma[given.ind, given.ind]

cMu <- drop(mu[req.ind] + C %*% solve(D) %*% (x.given - mu[given.ind]))

cVar <- B - C %*% solve(D) %*% t(C)

list(condMean=cMu, condVar=cVar)

}



n <- 10

A <- matrix(rnorm(n^2), n, n)

A <- A %*% t(A)

condNormal(x=c(1,1,0,0,-1), mu=rep(1,n), sigma=A, req=c(2,3,5),
given=c(1,4,7,9,10))



Best regards,

Ravi

Hi,

My previous version of the conditional MVN function had a bug in that it would
not work when conditional distribution was required for a single variable. I
fixed this and also made a few minor changes. Here is the new version.

condNormal <- function(x.given, mu, sigma, given.ind, req.ind){
# Returns conditional mean and variance of x[req.ind] 
# Given x[given.ind] = x.given
# where X is multivariate Normal with
# mean = mu and covariance = sigma
# 
B <- sigma[req.ind, req.ind]
C <- sigma[req.ind, given.ind, drop=FALSE]
D <- sigma[given.ind, given.ind]
CDinv <- C %*% solve(D)
cMu <- c(mu[req.ind] + CDinv %*% (x.given - mu[given.ind]))
cVar <- B - CDinv %*% t(C)
list(condMean=cMu, condVar=cVar)
}

n <- 10
A <- matrix(rnorm(n^2), n, n)
A <- A %*% t(A)
condNormal(x=c(1,1,0,0,-1), mu=rep(1,n), sigma=A, req=c(2,3,5),
given=c(1,4,7,9,10))
condNormal(x=c(1,1,0,0,-1), mu=rep(1,n), sigma=A, req=2, given=c(1,4,7,9,10))

As far as I know, there is nothing related to multivariate normal distributions
in "stats".  Hence, it seems like this function might be more useful
in a contributed package such as "fMultivar" or "mvtnorm".

Best,
Ravi.
________________________________________
From: r-devel-bounces at r-project.org [r-devel-bounces at r-project.org] on
behalf of Ravi Varadhan [rvaradhan at jhmi.edu]
Sent: Sunday, March 04, 2012 10:32 AM
To: r-devel at r-project.org
Subject: [Rd] Conditional means and variances of a multivariate normal
distribution

Hi,



Let X = (x_1, x_2, ... ,  x_p) be multivariate normal with mean, mu = (mu_1, ...
, mu_p) and covariance = Sigma.  I was looking for an R function to compute
conditional mean and conditional variance of a given subset of X given another
subset of X.  While this is trivially easy to do, there is nothing in
"base" for doing this, at least nothing that I am aware of.  I am also
not aware of anything in the contributed packages (although my search was not
comprehensive).  I feel that this would be a useful addition, if it is not
already there.  I have written this following function, which I am sure can be
improved a lot (including better argument names!). I would like to hear your
thought on this.



condNormal <- function(x.given, mu, sigma, req.ind, given.ind){

# Returns conditional mean and variance of x[req.ind]

# Given x[given.ind] = x.given

# where X is multivariate Normal with

# mean = mu and covariance = sigma

#

B <- sigma[req.ind, req.ind]

C <- sigma[req.ind, given.ind]

D <- sigma[given.ind, given.ind]

cMu <- drop(mu[req.ind] + C %*% solve(D) %*% (x.given - mu[given.ind]))

cVar <- B - C %*% solve(D) %*% t(C)

list(condMean=cMu, condVar=cVar)

}



n <- 10

A <- matrix(rnorm(n^2), n, n)

A <- A %*% t(A)

condNormal(x=c(1,1,0,0,-1), mu=rep(1,n), sigma=A, req=c(2,3,5),
given=c(1,4,7,9,10))



Best regards,

Ravi

______________________________________________
R-devel at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel