Inference methods comparison for the simplified Naomi model using maximum mean discrepancy

Abstract

Background We have run the simplified Naomi model using a range of inference methods: TMB, aghq and tmbstan.

Task In this report, we compare the accuracy of the posterior distributions obtained from these inference methods using maximum mean discrepancy.

1 Background

We compare the inference results from TMB, aghq, adam, and tmbstan. Import these inference results as follows:

tmb <- readRDS("depends/tmb.rds")
aghq <- readRDS("depends/aghq.rds")
tmbstan <- readRDS("depends/tmbstan.rds")

Check that the parameters (latent field, hyperparameters, model outputs) sampled from each of the four methods are the same:

stopifnot(names(tmb$fit$sample) == names(aghq$quad$sample))
stopifnot(names(tmb$fit$sample) == names(tmbstan$mcmc$sample))

2 Maximum mean discrepancy

par_samp_matrix <- function(sample) {
  x <- sample[!(stringr::str_ends(names(sample), "_out") | stringr::str_ends(names(sample), "_ll"))]
  do.call(rbind, x)
}

tmb_samples <- par_samp_matrix(tmb$fit$sample)
aghq_samples <- par_samp_matrix(aghq$quad$sample)
tmbstan_samples <- par_samp_matrix(tmbstan$mcmc$sample)

#' Keep only the final 1000 samples from tmbstan
tmbstan_samples <- tmbstan_samples[, 4001:5000]

mmd_tmb <- kernlab::kmmd(tmb_samples, tmbstan_samples, kernel = "rbfdot")

## Using automatic sigma estimation (sigest) for RBF or laplace kernel

mmd_aghq <- kernlab::kmmd(aghq_samples, tmbstan_samples, kernel = "rbfdot")

## Using automatic sigma estimation (sigest) for RBF or laplace kernel

#' Take sigma to be the mean estimated for the two different methods
kpar <- list(sigma = mean(mmd_tmb@kernelf@kpar$sigma, mmd_aghq@kernelf@kpar$sigma))

#' Rerun with sigma the same across the different tests
mmd_tmb <- kernlab::kmmd(tmb_samples, tmbstan_samples, kernel = "rbfdot", kpar = kpar, alpha = 0.05)
mmd_aghq <- kernlab::kmmd(aghq_samples, tmbstan_samples, kernel = "rbfdot", kpar = kpar, alpha = 0.05)

mmd_tmb

## Kernel Maximum Mean Discrepancy object of class "kmmd" 
##  
## Gaussian Radial Basis kernel function. 
##  Hyperparameter : sigma =  0.0012898572474155 
## 
## 
##  H0 Hypothesis rejected :  FALSE
##  Rademacher bound :  0.246480216773626
## 
##  1st and 3rd order MMD Statistics :  0.0801472436530942 0.00486463553339027

mmd_aghq

## Kernel Maximum Mean Discrepancy object of class "kmmd" 
##  
## Gaussian Radial Basis kernel function. 
##  Hyperparameter : sigma =  0.0012898572474155 
## 
## 
##  H0 Hypothesis rejected :  FALSE
##  Rademacher bound :  0.246480216773636
## 
##  1st and 3rd order MMD Statistics :  0.0713909472318657 0.00352025586168863

mmd_tmb@mmdstats

## [1] 0.080147244 0.004864636

mmd_aghq@mmdstats

## [1] 0.071390947 0.003520256

round(100 * (mmd_aghq@mmdstats[1] - mmd_tmb@mmdstats[1]) / mmd_tmb@mmdstats[1])

## [1] -11

round(100 * (mmd_aghq@mmdstats[2] - mmd_tmb@mmdstats[2]) / mmd_tmb@mmdstats[2])

## [1] -28