emix {bisque}

Compute expectations via weighted mixtures

Description

Approximates expectations of the form

E[h(\theta)] = \int h(\theta) f(\theta) d\theta

using a weighted mixture

E[h(\theta)] \approx \sum_{j=1}^k h(\theta^{(k)}) w_k

Usage

emix(h, params, wts, ncores = 1, errorNodesWts = NULL, ...)

Arguments

h

Function for which the expectation should be taken. The function should be defined so it is can be called via f(params, ...). Additional parameters may be passed to h via ....

params

Matrix in which each row contains parameters at which h should be evaluated. The number of rows in params should match the number of mixture components k.

wts

vector of weights for each mixture component

ncores

number of cores over which to evaluate mixture. this function assumes a parallel backend is already registered.

errorNodesWts

list with elements inds and weights that point out which params get used to compute an approximation of the quadrature error.

...

additional arguments to be passed to h

Examples

# density will be a mixture of betas
params = matrix(exp(2*runif(10)), ncol=2)

# mixture components are equally weighted
wts = rep(1/nrow(params), nrow(params))

# compute mean of distribution by cycling over each mixture component
h = function(p) { p[1] / sum(p) }

# compute mixture mean
mean.mix = emix(h, params, wts)

# (comparison) Monte Carlo estimate of mixture mean
nsamples = 1e4
component = sample(x = 1:length(wts), size = nsamples, prob = wts, 
                   replace = TRUE)
x = sapply(component, function(cmp) {
  rbeta(n = 1, shape1 = params[cmp, 1], shape2 = params[cmp, 2])
})
mean.mix.mc = mean(x)

# compare estimates
c(emix = mean.mix, MC = mean.mix.mc)