############################################################################
#         R FUNCTIONS FOR THE MULTIVARIATE NORMAL DISTRIBUTION
#
#         (by Jeffrey S. Rosenthal, probability.ca, 2009)
# 
# The base package of the statistical software R does not seem to include
# built-in functions for the multivariate normal distribution.  So, I provide
# them here.  (Released under GPL.)
############################################################################

# dmvn: density function of the multivariate normal distribution
dmvn = function(x, mean = rep(0,length(x)), cov = diag(length(mean)),
							log = FALSE) {
  thedim = length(mean);
  if (!is.matrix(cov))
    stop("Error: covariance argument must be a matrix.");
  if ( (!is.real(cov)) || (!isSymmetric(cov)) )
    stop("Error: covariance matrix must be real and symmetric.");
  if ( (nrow(cov) != thedim) || (nrow(cov) != thedim) )
    stop("Error: mean and covariance matrix have incompatible sizes.");
  if (length(x) != thedim)
    stop("Error: x and mean have different lengths.");
  logresult = (-thedim/2) * log(2*3.1415926) - 0.5 * log(abs(det(cov))) -
		0.5 * t(x-mean) %*% solve(cov) %*% (x-mean) ;
  logresult = logresult[1,1];
  if (log)
    return(logresult)
  else
    return(exp(logresult));
}


# rmvn: sample a vector of multivariate normal random variables
# (Note: at least one of mean and cov must be specified.)
rmvn = function(mean = rep(0,nrow(cov)), cov = diag(length(mean)) ) {
  thedim = length(mean);
  if (!is.matrix(cov))
    stop("Error: covariance argument must be a matrix.");
  if ( (!is.real(cov)) || (!isSymmetric(cov)) )
    stop("Error: covariance matrix must be real and symmetric.");
  if ( (nrow(cov) != thedim) || (nrow(cov) != thedim) )
    stop("Error: mean and covariance matrix have incompatible sizes.");
  A = t(chol(cov));   # this ensures that: A A^T = cov
  return( t(mean + A %*% rnorm(thedim))[1,] );
}


# Example:
mycov = matrix( c(1,0.98,0.98,1), nrow=2 )
mymean = c(10,100)
x = c(11, 101)
print( rmvn(mymean,mycov) )
print( dmvn(x,mymean,mycov) )