# 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) } }