Generate a pseudo-Hessian matrix from MCMC series
The Hessian matrix at a point P (P is a vector) is the matrix of partial derivatives at this point. It is often used after non-linear optimization and the standard error of parameters can be obtained using the square-root of the diagonal of the inverse of the Hessian matrix. But when you are estimate SE like this, you lost the covariances between parameters. Then if you estimate random numbers using the SE of parameters, the original structure is lost. It is much better to estimate random numbers directly from the Hessian matrix. Then you keep the variance-covariance structure. When the model is fitted using MCMC, no Hessian matrix is directly available. A pseudo-Hessian matrix can be calculated using: hessian <- solve(cov( [matrix of values for different iterations] )) Let do an example with the fit of a normal distribution: First using maximum likelihood: val <- rnorm(30, 10, 2) library(MASS) ft <- fitdistr(val, densfun="normal") The Hessian matrix is not available