In sieste/doit: Bayesian Computation using Design of Experiments-Based Interpolation Technique

knitr::opts_chunk$set(
  fig.path = paste('figure/', knitr::current_input(), '-figs/', sep=''),
  fig.width = 7,
  fig.height = 5
)

library(doit)
library(lhs)
library(tidyverse)

suppressPackageStartupMessages(library(doit))
suppressPackageStartupMessages(library(lhs))
suppressPackageStartupMessages(library(tidyverse))

The target distribution

We use the example model from Joseph (2012), namely a Bernoulli likelihood with logit link function and a Normal prior:

• $r \sim N(\mu, \tau^2)$
• $p = \frac{1}{1+exp(-r)}$
• $y \sim Ber(p)$
• where $y=1$, $\mu = 1$, $\tau = 4$

The unnormalised posterior is given by

h = function(r) {
  p = 1/(1+exp(-r))
  dbinom(1, 1, p) * dnorm(r, 1, 4)
}

Design points

We evaluate the unnormalised posterior at 10 equally spaced design points:

design = data.frame(x=seq(-10, 20, len=10))
design$f = sapply(design$x, h)

Estimate the optimal kernel width

The DoIt method works by approximating the target density by a weighted sum of Normal distributions. The optimal variance parameter w of the Normal distributions is estimated by the function doit_estimate_w():

w = doit_estimate_w(design)
print(w)

Evaluate the DoIt approximation

The function doit_fit() estimates the parameters of the DoIt method and saves them in an object of class doit:

doit = doit_fit(design, w=w)

We specify the parameter values at which we want to approximate the target posterior and use doit_approx() to calculate the approximation:

theta_eval = data.frame(x = seq(-15,25,.1))
doit_out = doit_approx(doit, theta_eval)

The approximation function returns the approximated density value and the estimation error variance:

head(doit_out)

The output of doit_approx() works nicely with the ggplot2 package:

ggplot() + 
geom_line(data=doit_out, aes(x=x, y=dens_approx, colour=variance)) +
geom_point(data=design, aes(x=x, y=0), pch=1, cex=2)

Add more design points

To extend the existing design, we use doit_propose_new() to choose a new design point where the current estimation variance is particularly high:

x_n = doit_propose_new(doit)
design_new = mutate(x_n, f=h(x))

Instead of using doit_fit() again on the original design plus additional new point, we can use the doit_update() function which is more efficient. Note that the kernel width w is not updated:

doit = doit_update(doit, design_new)

We plot the new approximation and highlight the added design point:

doit_out = doit_approx(doit, theta_eval)
ggplot() + 
geom_point(data=design_new, aes(x=x, y=0), pch=1, cex=2, col='red') +
geom_point(data=design, aes(x=x, y=0), pch=1, cex=2) +
geom_line(data=doit_out, aes(x=x, y=dens_approx, colour=variance))

Posterior expectation and variance

We have functions to approximate posterior mean and variance:

doit_expectation(doit)
doit_variance(doit)

sieste/doit documentation built on May 9, 2019, 4:10 p.m.