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 in ft list, but the variance-covariance matrix is available:

ft$vcov

          mean         sd

mean 0.1033382 0.00000000

sd   0.0000000 0.05166912

Then it is easy to get the Hessian matrix:

solve(ft$vcov)

         mean       sd

mean 9.676959  0.00000

sd   0.000000 19.35392

To be sure, let estimate the true Hessian matrix:

library("numDeriv")

dnormx <- function(x, data) {

  data <- unlist(data)

  return(-sum(dnorm(data, mean=x['mean'], sd=x['sd'], log=TRUE)))

}

hessian(dnormx, x=ft$estimate, data=val)

              [,1]          [,2]

[1,]  9.676959e+00 -6.274874e-10

[2,] -6.274874e-10  1.935392e+01

Perfect. Now let try with MCMC:

library(HelpersMG)
require(coda)
parameters_mcmc <- data.frame(Density=c('dnorm', 'dlnorm'), 
                              Prior1=c(10, 0.5), Prior2=c(2, 0.5), SDProp=c(0.35, 0.2), 
                              Min=c(-3, 0), Max=c(100, 10), Init=c(10, 2), stringsAsFactors = FALSE, 
                              row.names=c('mean', 'sd'))
# Use of trace and traceML parameters
# trace=1 : Only one likelihood is printed
mcmc_run <- MHalgoGen(n.iter=50000, parameters=parameters_mcmc, data=val, 
                      likelihood=dnormx, n.chains=1, n.adapt=100, thin=1, trace=1)

The variance-covariance matrix can be calculated easily:

cov(mcmc_run$resultMCMC[[1]])
              mean            sd
mean  0.1100181574 -0.0007734971
sd   -0.0007734971  0.0571208557

And the pseudo-Hessian is then:

solve(cov(mcmc_run$resultMCMC[[1]]))
          mean         sd
mean 9.0902742  0.1230952
sd   0.1230952 17.5084074

This pseudo-Hessian can be used to calculate random numbers using the function RandomFromHessianOrMCMC() from HelpersMG package:

r <- RandomFromHessianOrMCMC(Hessian = solve(cov(mcmc_run$resultMCMC[[1]])),   
fitted.parameters = c(mean=mean(mcmc_run$resultMCMC[[1]][, "mean"]), 
sd=mean(mcmc_run$resultMCMC[[1]][, "sd"])))
Estimation using variance-covariance matrix

And if we look at the 10000 random numbers, they respect well the averages:

colMeans(r$random)

     mean        sd 

10.226867  1.824722 

But also the covariance structure:

cov(as.matrix(r$random))

              mean            sd

mean  0.1113069709 -0.0007618878

sd   -0.0007618878  0.0585928419

And of course the Hessian matrix:

solve(cov(as.matrix(r$random)))

          mean         sd

mean 8.9849630  0.1168323

sd   0.1168323 17.0684504