gompertz: Posterior distribution (and related quantities) for the...

In erlisR/iLaplaceExamples: Extra Utilities and Examples for iLaplace

Description

Function nlpost_gomp gives the non-normalised negative log-posterior distribution for the Gompertz distribution, grad_gomp gives its gradient and hess_gomp gives the Hessian matrix. The parameter of interest is (exp(α), exp(β)), with (α,β) being reals. See Details for more information.

Usage

nlpost_gomp(param, data)

grad_gomp(param, data)

hess_gomp(param, data)

Arguments

param	The bivariate vector of parameters.
data	The vector of data.

Details

The prior distribution is bivariate normal with independent components, with mean (0,0) and variance (100, 100).

Value

double

bivariate vector

2 by 2 matrix

See Also

MHmcmc

Examples

## Not run: 

library(VGAM)
library(iLaplace)
library(iLaplaceExamples)
set.seed(12)

# generate the data
data <- rgompertz(n = 50, shape = exp(2), scale = exp(3))

ff <- function(x, ...) nlpost_gomp(param = x, ...)
ff.gr <- function(x, ...) grad_gomp(param = x, ...)
ff.hess <- function(x, ...) hess_gomp(param = x, ...)

# find the posterior mode and Hessian at the mode
opt.post = nlminb(c(2,3), obj = ff, gradient = ff.gr,
                  hessian = ff.hess, data = data)
opt.post$hessian = ff.hess(opt.post$par, data = data)

# draw a posterior sample
mcmc.gomp = MHmcmc(logfun = function(x) -ff(x, data),
                   burnin = 5000, mcmc = 1e+6, thin = 2, tune = 1.1,
                   V = solve(opt.post$hessian), df = 5,
                   theta.init = opt.post$par, verbose = 10000)
# look at the plots
plot(mcmc.gomp)

# marginal likelihood by the improved Laplace of Ruli et al. (2015)
iLap.logM <- iLaplace(fullOpt = opt.post, ff = ff, ff.gr = ff.gr, ff.hess = ff.hess,
                      control = list(sp.points = 100, delta = 13, n.cores = 2),
                      clEvalQ = c("iLaplaceExamples"), data = data)

# marginal likelihood by the standard Laplace
Lap.logM <- log(2*pi) - opt.post$objective - 0.5*determinant(opt.post$hessian)$mod

# marginal likelihood by improtance sampling
is.logM <- isML(logfun = function(x) -ff(x, data), nsim = 1e+6,
                theta.hat = opt.post$par, tune = 1.2, V = solve(opt.post$hessian),
               df = 5, verbose = 100000)

## End(Not run)

erlisR/iLaplaceExamples documentation built on May 16, 2019, 8:48 a.m.