wMix {bisque} | R Documentation |
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.
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
|
f2 |
evaluates |
w |
|
f1.precompute |
function that pre-computes parameters for evaluating
|
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 |
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(\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)
.
# 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)))