calcBF.R

# I assume that you have already run MCMC, and parameter vectors
# are stored in a n-by-p matrix named 'mcmcParameterEstimates'.  The
# posterior modes of each parameter (the mode of each column) are
# stored in a vector named 'parModes'.


# calculate the covariance matrix from parameter estimates
covMat <- cov(mcmcParameterEstimates)

# scale the covariance matrix estimated from mcmc samples (optional)
scalar <- 1.0;
covMat <- covMat * scalar;
# use the 'try' wrapper to not show warnings and invert the matrix
INV <- try(solve(covMat), silent=TRUE)

# logW stores the logarithm of the importance weights   
logW <- 0
if(is.null(dim(INV)) == TRUE) { 
   # no solution (matrix may be singular), set everything equal to NA
   logW[1:numMCMCSamples] <- NA;
}
if(is.null(dim(INV)) != TRUE) {
   # calculate the logarithm of the determinant
   LOGDET <- sum(logb(abs(diag(qr(covMat)[[1]]))))
   # cycle through each parameter vector, calculating the importance weight
   for(i in 1:numMCMCSamples) {
      logW[i] <- fastdmnorm(mcmcParameterEstimates[i,], mean=as.vector(parModes), INV, LOGDET)
   }
}