Description Usage Arguments Details Value See Also Examples
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.
1 2 3 4 5 |
param |
The bivariate vector of parameters. |
data |
The vector of data. |
The prior distribution is bivariate normal with independent components, with mean (0,0) and variance (100, 100).
double
bivariate vector
2 by 2 matrix
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 | ## 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)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.