Introduction

Let X be a vector of parameters from a MVN(X, V) draw, the importance ratio is the following:

\[\frac{L(X | data)}{\text{MVN}( X | \bf{X}, \bf{V})}\]

Gaussian example

First simulation data

n <- 1e3
set.seed(1129)
dat <- rnorm(n)

Fit the data via mle

library(bbmle)
## Loading required package: stats4
nll <- function(m=0,s=1){
  return(-sum(dnorm(dat,m,s,log=TRUE)))
}
  
mlefit <- mle2(nll)

print(coef(mlefit))
##          m          s 
## 0.01242379 0.97781068
print(vcov(mlefit))
##              m            s
## m 9.561137e-04 5.323910e-11
## s 5.323910e-11 4.780576e-04

We have two methods for the covariance matrix. The book suggest the pseudo-variance using the generalized cholesky of the generalized inverse of the hessian.

samp_size <- 1e3

hh <- mlefit@details$hessian
V <- bdsmatrix:::gchol(MASS:::ginv(hh))
tV <- matrix(V@.Data,nrow=2)
V <- t(tV)

vv <- tV %*% V
mv_samps0 <- MASS:::mvrnorm(samp_size, mu = coef(mlefit), Sigma=vv)

print(summary(mv_samps0))
##        m                  s         
##  Min.   :0.009425   Min.   :0.9765  
##  1st Qu.:0.011825   1st Qu.:0.9775  
##  Median :0.012424   Median :0.9778  
##  Mean   :0.012433   Mean   :0.9778  
##  3rd Qu.:0.013036   3rd Qu.:0.9781  
##  Max.   :0.015726   Max.   :0.9791
vv <- vcov(mlefit)
mv_samps <- MASS:::mvrnorm(samp_size, mu = coef(mlefit), Sigma=vv)

print(summary(mv_samps))
##        m                   s         
##  Min.   :-0.083534   Min.   :0.8959  
##  1st Qu.:-0.009027   1st Qu.:0.9635  
##  Median : 0.010543   Median :0.9784  
##  Mean   : 0.010896   Mean   :0.9786  
##  3rd Qu.: 0.033203   3rd Qu.:0.9931  
##  Max.   : 0.140189   Max.   :1.0536

We are just going to stick with the estimated covariance matrix.

Posterior_negLL <- sapply(1:samp_size,function(x){nll(mv_samps[x,1],mv_samps[x,2])})
approx_LL <- sapply(1:samp_size
  , function(x){
    mvtnorm::dmvnorm(mv_samps[x,]
      , mean = coef(mlefit)
      , sigma=vv
      , log = TRUE
      )
    }
)

Log_imp_wts <- -Posterior_negLL - approx_LL
print(summary(Log_imp_wts))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   -1403   -1402   -1402   -1402   -1402   -1401
Log_scaled_imp_wts <- Log_imp_wts - max(Log_imp_wts)
print(summary(Log_scaled_imp_wts))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -2.0226 -0.8352 -0.8200 -0.8209 -0.7987  0.0000
imp_wts <- exp(Log_scaled_imp_wts)
imp_wts <- imp_wts/sum(imp_wts) 
print(summary(imp_wts))
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0002990 0.0009804 0.0009955 0.0010000 0.0010169 0.0022602

Diagnostics

We can look at pairs plot, effective sample size and weighted quantiles.

pp <- cbind(mv_samps,-Posterior_negLL, approx_LL)
print(pairs(pp))

## NULL
eff_samp <- 1/sum(imp_wts^2)
print(eff_samp) 
## [1] 989.2847

Looks good!

wq <- sapply(1:2
  , function(x){Hmisc::wtd.quantile(mv_samps[,x]
    , weights = imp_wts
    , probs = c(0.025, 0.975)
    , normwt = TRUE
    )
  }
)
print(t(wq))
##             2.5%     97.5%
## [1,] -0.05107104 0.0710281
## [2,]  0.93863089 1.0269053
print(confint(mlefit))
##         2.5 %     97.5 %
## m -0.04823918 0.07308665
## s  0.93647673 1.02228278

Close enough.