wMix {bisque}

Construct a weighted mixture object

Description

For a Bayesian model

X ~ f(X | \theta_1, \theta_2)

(\theta_1, \theta_2) ~ f(\theta_1, \theta_2),

the marginal posterior f(\theta_1 | X) distribution can be approximated via weighted mixtures via

f(\theta_1 | X) \approx \sum_{j=1}^K f(\theta_1 | X, \theta_2) w_j

where w_j is based on f(\theta_2^{(j)} | X) and weights \tilde w_j, where \theta_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(\theta_2|X), after transforming \theta_2 to an unconstrained support. For best results, \theta_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(\theta_1 | X, \theta_2). f1 must be able to be called via f1(theta1, params, log, ...).

theta1

a matrix of parameters at which to evaluate f(\theta_1 | X, \theta_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(\theta_1 | X, \theta_2). In most cases params will equal theta_2, but in some cases, f(\theta_1 | X, \theta_2) depends on functions of \theta_2, which can be pre-evaluated as the weighted mixture approximation is constructed.

log

TRUE to return ln(f(\theta_1 | X, \theta_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(\theta_2| X) and wrapper functions to generate quadrature grid.

f1.precompute

function that pre-computes parameters for evaluating f(\theta_1 | X, \theta_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(\theta_1|\theta_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(\theta_1 | \theta_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(\theta_1|X). f is callable via f(theta1, log, ...).

mix

A matrix containing the pre-computed parameters for evaluating the mixture components f(\theta_1 | \theta_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(\theta_2|X). However, posterior expectations of f(\theta_1|X) can also be computed when expectations of f(\theta_1|\theta_2, X) are known. The elements of expectation are

Eh

Function to compute E[h(\theta_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 \theta_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 \theta_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(\theta_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)))