wMix {bisque}  R Documentation 
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 sparsegrid quadrature integration scheme. The quadrature rule is developed by finding the posterior mode of f(θ_2X), 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.
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, ... )
f1 
evaluates f(θ_1  X, θ_2).

f2 
evaluates f(theta_2  X). 
w 

f1.precompute 
function that precomputes parameters for evaluating
f(θ_1  X, θ_2). 
spec 
Specification of whether 
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 
c.level 
accuracy level of the numerical approximation for 
c.init 
initial guess for mode of 
c.link 
character vector that specifies transformations used during
optimization and integration of 
c.link.params 
Optional list of additional parameters for link
functions used with 
c.optim.control 
Arguments used to find mode of 
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 
... 
Additional arguments to pass to 
A list with class wMix
, which contains the following items.
f
Function for evaluating the posterior density
f(θ_1X). f
is callable via
f(theta1, log, ...)
.
mix
A matrix containing the precomputed 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(θ_2X). However, posterior
expectations of f(θ_1X) 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 sparsequadrature 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(θ_1X).
# Use BISQuE to approximate the marginal posterior distribution for unknown # population f(Nc, r) for the fur seals capturerecapture 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(1alpha) 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(1alpha[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 = Nr, 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) ((1p)*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)))