wBuild: Derive parameters for building integration grids

In jmhewitt/bisque: Approximate Bayesian Inference via Sparse Grid Quadrature Evaluation (BISQuE) for Hierarchical Models

Description

Note: w is defined on the transformed scale, but for convenience f is defined on the original scale.

Usage

wBuild(
  f,
  init,
  dim.theta2 = length(init),
  approx = "gaussian",
  link = rep("identity", length(init)),
  link.params = rep(list(NA), length(init)),
  optim.control = list(maxit = 5000, method = "BFGS"),
  ...
)

Arguments

f	function used to derive the weight function w. f must be able to be called via f(par, log, ...)
init	initial guess for mode of f.
dim.theta2	wBuild assumes par is partitioned such that par=c(theta1,theta2). dim.theta2 specifies the size of the partition. The default is to assume that f is defined without a theta1 component.
approx	Style of approximation (i.e., w) to be created from mode of f. 'gaussian' Gaussian approximation for theta2 at the mode of f. Assumes f is proportional to the marginal posterior density for theta2. 'condgauss' Gaussian approximation for theta2 at the mode of f. The approximation is conditioned on the value of the mode for theta1. Assumes f is proportional to the joint posterior density for theta1,theta2. 'condgauss-laplace' Gaussian approximation for theta2 at the mode of f. The approximation is conditioned on a separate laplace approximation of the marginal posterior mode for theta1. Assumes f is proportional to the joint posterior density for theta1,theta2. 'margauss' Gaussian approximation for theta2 at the mode of f. Assumes f is proportional to the joint posterior density for theta1,theta2., then uses the marginal mean and covariance from the posterior's gaussian approximation.
link	character vector that specifies transformations used during optimization and integration of f(θ_2 | X). While θ_2 may be defined on arbitrary support, wtdMix performs optimization and integration of θ_2 on an unconstrained support. The link vector describes the transformations that must be applied to each element of θ_2. Jacobian functions for the transformations will automatically be added to the optimization and integration routines. Currently supported link functions are 'log', 'logit', and 'identity'.
link.params	Optional list of additional parameters for link functions. For example, the logit function can be extended to allow mappings to any closed interval. There should be one list entry for each link function. Specify NA if no additional arguments are passed.
optim.control	List of arguments to pass to stat::optim when used to find mode of f. maxit Maximum number of iterations to run optim for. method Optimization routine to use with optim.
...	additional arguments needed for function evaluation.

Examples

# Use BISQuE to approximate the marginal posterior distribution for unknown
# population f(N|c, r) for the fur seals capture-recapture data example in 
# Givens and Hoeting (2013), example 7.10.

data('furseals')

# define theta transformation and jacobian
tx.theta = function(theta) { 
  c(log(theta[1]/theta[2]), log(sum(theta[1:2]))) 
}
itx.theta = function(u) { 
  c(exp(sum(u[1:2])), exp(u[2])) / (1 + exp(u[1])) 
}
lJ.tx.theta = function(u) {
  log(exp(u[1] + 2*u[2]) + exp(2*sum(u[1:2]))) - 3 * log(1 + exp(u[1]))
}

# compute constants
r = sum(furseals$m)
nC = nrow(furseals)

# set basic initialization for parameters
init = list(U = c(-.7, 5.5))
init = c(init, list(
  alpha = rep(.5, nC),
  theta = itx.theta(init$U),
  N = r + 1
))


post.alpha_theta = function(theta2, log = TRUE, ...) {
  # Function proportional to f(alpha, U1, U2 | c, r) 
  
  alpha = theta2[1:nC]
  u = theta2[-(1:nC)]
  theta = itx.theta(u)
  p = 1 - prod(1-alpha)
  
  res = - sum(theta)/1e3 - r * log(p) + lJ.tx.theta(u) - 
    nC * lbeta(theta[1], theta[2])
  for(i in 1:nC) {
    res = res + (theta[1] + furseals$c[i] - 1)*log(alpha[i]) + 
      (theta[2] + r - furseals$c[i] - 1)*log(1-alpha[i])
  }
  
  if(log) { res } else { exp(res) }
}

post.N.mixtures = function(N, params, log = TRUE, ...) {
  # The mixture component of the weighted mixtures for f(N | c, r)
  dnbinom(x = N-r, size = r, prob = params, log = log)
}

mixparams.N = function(theta2, ...) {
  # compute parameters for post.N.mixtures
  1 - prod(1 - theta2[1:nC])
}


w.N = wBuild(f = post.alpha_theta, init = c(init$alpha, init$U), 
             approx = 'gauss', link = c(rep('logit', nC), rep('identity', 2)))

m.N = wMix(f1 = post.N.mixtures, f1.precompute = mixparams.N, 
           f2 = post.alpha_theta, w = w.N)



# compute posterior mean
m.N$expectation$Eh.precompute(h = function(p) ((1-p)*r/p + r), 
                                   quadError = TRUE)

# compute posterior density
post.N.dens = data.frame(N = r:105)
post.N.dens$d = m.N$f(post.N.dens$N)

# plot posterior density
plot(d~N, post.N.dens, ylab = expression(f(N~'|'~bold(c),r)))

jmhewitt/bisque documentation built on Feb. 9, 2020, 2:36 a.m.