# 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

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 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, ...). theta1a 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 paramsa 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. logTRUE 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

wts
  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 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 # 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)))