wMix: Construct a weighted mixture object

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

Description

For a Bayesian model

X ~ f(X | θ_1, θ_2)

(θ_1, θ_2) ~ f(θ_1, θ_2),

the marginal posterior f(θ_1 | X) distribution can be approximated via weighted mixtures via

f(θ_1 | X) \approx ∑_{j=1}^K f(θ_1 | X, θ_2) w_j

where w_j is based on f(θ_2^{(j)} | X) and weights \tilde w_j, where θ_2^{(j)} and \tilde w_j are nodes and weights for a sparse-grid quadrature integration scheme. The quadrature rule is developed by finding the posterior mode of f(θ_2|X), after transforming θ_2 to an unconstrained support. For best results, θ_2 should be a continuous random variable, or be able to be approximated by one.

Usage

wMix(
  f1,
  f2,
  w,
  f1.precompute = function(x, ...) {     x },
  spec = "ff",
  level = 2,
  c.int = NULL,
  c.level = 2,
  c.init = NULL,
  c.link = rep("identity", length(c.init)),
  c.link.params = rep(list(NA), length(c.init)),
  c.optim.control = list(maxit = 5000, method = "BFGS"),
  ncores = 1,
  quadError = TRUE,
  ...
)

Arguments

f1	evaluates f(θ_1 | X, θ_2). f1 must be able to be called via f1(theta1, params, log, ...). theta1 a matrix of parameters at which to evaluate f(θ_1 | X, θ_2). each row should be one set of values at which the density should be evaluated params a vector of parameters needed to evaluate f(θ_1 | X, θ_2). In most cases params will equal theta_2, but in some cases, f(θ_1 | X, θ_2) depends on functions of θ_2, which can be pre-evaluated as the weighted mixture approximation is constructed. log TRUE to return ln(f(θ_1 | X, θ_2)) ... additional arguments needed for function evaluation
f2	evaluates f(theta_2 | X). f2 must be able to be called via f2(theta2, log, ...).
w	wBuild object created by wBuild function. w contains posterior mode of f(θ_2| X) and wrapper functions to generate quadrature grid.
f1.precompute	function that pre-computes parameters for evaluating f(θ_1 | X, θ_2). f1.precompute must be able to be called via f1.precompute(theta2, ...) and return the argument params for the function f1.
spec	Specification of whether f1 and f2 are known exactly, or need numerical approximation to determine integration constants. 'ff' if both functions are known, 'gg' if f1 is proportional to the full conditional distribution f(θ_1|θ_2,X), but needs the integration constant computed, and if the marginal posterior f(theta_2|X) is equal to f2 times the integration constant that needs to be numerically approximated.
level	accuracy level of the numerical approximation (typically number of grid points for the underlying 1D quadrature rule) [description from mvQuad::createNIGrid]
c.int	If spec=='gg', then c.int specifies the function that can be integrated in order to yield the missing integration constant.
c.level	accuracy level of the numerical approximation for c.int (typically number of grid points for the underlying 1D quadrature rule) [description from mvQuad::createNIGrid]
c.init	initial guess for mode of c.int.
c.link	character vector that specifies transformations used during optimization and integration of c.int. See corresponding documentation in wBuild function for more details.
c.link.params	Optional list of additional parameters for link functions used with c.int. See corresponding documentation in wBuild function for more details.
c.optim.control	Arguments used to find mode of c.int. See corresponding documentation in wBuild function for more details.
ncores	number of cores used to parallelize computation of parameters for f(θ_1 | θ_2, X).
quadError	TRUE if integration nodes and weight should be computed for the level-1 integration grid, so that quadrature approximation error can be estimated.
...	Additional arguments to pass to f1, f1.precompute, f12, and f2.

Value

A list with class wMix, which contains the following items.

f

Function for evaluating the posterior density f(θ_1|X). f is callable via f(theta1, log, ...).

mix

A matrix containing the pre-computed parameters for evaluating the mixture components f(θ_1 | θ_2, X). Each row of the matrix contains parameters for one of the K mixture components.

wts

Integration weights for each of the mixture components. Some of the weights may be negative.

expectation

List containing additional tools for computing posterior expectations of f(θ_2|X). However, posterior expectations of f(θ_1|X) can also be computed when expectations of f(θ_1|θ_2, X) are known. The elements of expectation are

Eh

Function to compute E[h(θ_2)|X]. Eh is callable via Eh(h, ...), where h is a function callable via h(theta2, ...) and ... are additional arguments to the function. The function h is evaluated at the quadrature nodes θ_2^{(j)}.

Eh.precompute

Exactly the same idea as Eh, but the function h is evalauted at the quadrature nodes after being passed through the function f1.precompute.

grid

The sparse-quadrature integration grid used. Helpful for seeing the quadrature nodes θ_2^{(j)}.

wts

The integration weights for approximating the expectation E[h]. Note that these integration weights may differ from the main integration weights for evaluating the posterior density f(θ_1|X).

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.