Create a centered and scaled sparse integration grid

Description

Enhances mvQuad::createNIGrid by shifting and scaling a sparse integration grid, and evaluating the weight function at each of the grid nodes.

Usage

createLocScaleGrid(
  mu = 0,
  prec = 1,
  level = 2,
  quadError = FALSE,
  prec.chol = chol(prec)
)

Arguments

See Also

mvQuad::createNIGrid

Examples

g = createLocScaleGrid(mu = c(1,0), prec = diag(c(1,.5)), level = 2 )

Evaluate a mixture density

Description

Evaluates mixture densities of the form

f(x) = \sum_{j=1}^k f(x|\theta^{(k)}) w_k

where the w_k are (possibly negative) weights that sum to 1 and f(x|\theta^{(k)}) are densities that are specified via parameters \theta^{(k)}, which are passed in the function argument params. A unique feature of this function is that it is able to evaluate mixture densities in which some of the mixture weights w_k are negative.

Usage

dmix(x, f, params, wts, log = FALSE, errorNodesWts = NULL, ...)

Arguments

Examples

# evaluate mixture density at these locations
x = seq(0, 1, length.out = 100)

# density will be a mixture of beta distributions
f = function(x, theta, log = FALSE) {
  dbeta(x, shape1 = theta[1], shape2 = theta[2], log = log)
}

# beta parameters are randomly assigned
params = matrix(exp(2*runif(10)), ncol=2)

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

# evaluate mixture density
fmix = dmix(x = x, f = f, params = params, wts = wts)

# plot mixture density
plot(x, fmix, type='l', ylab = expression(f(x)), 
     ylim = c(0, 4))

# plot component densities
for(i in 1:length(wts)){
  curve(f(x, params[i,]), col = 2, add = TRUE)
}

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

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)

Data from a capture-recapture study of fur seal pups

Description

These data are used in the book "Computational Statistics" by G.H. Givens and J.A. Hoeting (2013). They are discussed in Chapter 7, Examples 7.2,7.3,7.8, and Exercise 7.2.

Usage

data(furseals)

Format

A data.frame with variables:

i

The census attempt

c

Number of pups captured in census attempt

m

Number of newly captured pups

Details

As described by the authors:

Source: Richard Barker, University of Otago, New Zealand

Description: Data from a capture-recapture study conducted on the Otago Penninsula, South Island, New Zealand. Fur seal pups were marked and released during 7 census attempts in one season. The population is assumed closed. For each census attempt, the number of pups captured and the number of these captures corresponding to pups never previously caught are recorded.

Source

https://www.stat.colostate.edu/computationalstatistics/

https://www.stat.colostate.edu/computationalstatistics/datasets.zip

Examples

data("furseals")
str(furseals)

Named inverse transformation functions

Description

Evaluates the inverse of the named link function at the locations x.

Usage

itx(x, link, linkparams)

Arguments

Examples

bisque:::itx(0, 'logit', list(NA))

Jacobian for exponential transform

Description

Let X=exp(Y) be a transformation of a random variable Y. This function computes the jacobian J(x) when using the density of Y to evaluate the density of X via

f(x) = f_y(ln(x)) J(x)

where

J(x) = d/dx ln(x).

Usage

jac.exp(x, log = TRUE)

Arguments

Examples

jac.exp(1)

Jacobian for logit transform

Description

Let X=logit^{-1}(Y) be a transformation of a random variable Y. This function computes the jacobian J(x) when using the density of Y to evaluate the density of X via

f(x) = f_y(logit(x)) J(x)

where

J(x) = d/dx logit(x).

Usage

jac.invlogit(x, log = TRUE)

Arguments

Examples

jac.invlogit(1)

Jacobian for log transform

Description

Let X=log(Y) be a transformation of a random variable Y. This function computes the jacobian J(x) when using the density of Y to evaluate the density of X via

f(x) = f_y(exp(x)) J(x)

where

J(x) = d/dx exp(x).

Usage

jac.log(x, log = TRUE)

Arguments

Examples

jac.log(1)

Jacobian for logit transform

Description

Let X=logit(Y) be a transformation of a random variable Y that lies in the closed interval (L,U). This function computes the jacobian J(x) when using the density of Y to evaluate the density of X via

f(x) = f_y(logit^{-1}(x) * (U-L) + L) J(x)

where

J(x) = (U-L) d/dx logit^{-1}(x).

Usage

jac.logit(x, log = TRUE, range = c(0, 1))

Arguments

Examples

jac.logit(1)

Use sparse grid quadrature techniques to integrate (unnormalized) densities

Description

This function integrates (unnormalized) densities and may be used to compute integration constants for unnormalized densities, or to marginalize a joint density, for example.

Usage

kCompute(
  f,
  init,
  method = "BFGS",
  maxit = 10000,
  level = 2,
  log = FALSE,
  link = NULL,
  linkparams = NULL,
  quadError = FALSE,
  ...
)

Arguments

Examples

kCompute(dgamma, init = 1, shape=2, link='log', level = 5)

Wrapper to abstractly evaluate log-Jacobian functions for transforms

Description

Wrapper to abstractly evaluate log-Jacobian functions for transforms

Usage

logjac(x, link, linkparams)

Arguments

See Also

jac.log, jac.logit

Examples

bisque:::logjac(1, 'logit', list(NA))

Merge pre-computed components of f(theta1 | theta2, X)

Description

For use in the parallel call in wtdMix()

Usage

mergePars(x, y)

Arguments

Fit a spatially mean-zero spatial Gaussian process model

Description

Uses a Gibbs sampler to estimate the parameters of a Matern covariance function used to model observations from a Gaussian process with mean 0.

Usage

sFit(
  x,
  coords,
  nSamples,
  thin = 1,
  rw.initsd = 0.1,
  inits = list(),
  C = 1,
  alpha = 0.44,
  priors = list(sigmasq = list(a = 2, b = 1), rho = list(L = 0, U = 1), nu = list(L = 0,
    U = 1))
)

Arguments

Examples

library(fields)

simulate.field = function(n = 100, range = .3, smoothness = .5, phi = 1){
  # Simulates a mean-zero spatial field on the unit square
  #
  # Parameters:
  #  n - number of spatial locations
  #  range, smoothness, phi - parameters for Matern covariance function
  
  coords = matrix(runif(2*n), ncol=2)
  
  Sigma = Matern(d = as.matrix(dist(coords)), 
                 range = range, smoothness = smoothness, phi = phi)
  
  list(coords = coords,
       params = list(n=n, range=range, smoothness=smoothness, phi=phi),
       x = t(chol(Sigma)) %*%  rnorm(n))
}

# simulate data
x = simulate.field()

# configure gibbs sampler  
it = 100

# run sampler using default posteriors
post.samples = sFit(x = x$x, coords = x$coords, nSamples = it)

# build kriging grid
cseq = seq(0, 1, length.out = 10)
coords.krig = expand.grid(x = cseq, y = cseq)

# sample from posterior predictive distribution
burn = 75
samples.krig = sKrig(x$x, post.samples, coords.krig = coords.krig, burn = burn)

Draw posterior predictive samples from a spatial Gaussian process model

Description

Draw posterior predictive samples from a spatial Gaussian process model

Usage

sKrig(x, sFit, coords.krig, coords = sFit$coords, burn = 0, ncores = 1)

Arguments

Examples

library(fields)

simulate.field = function(n = 100, range = .3, smoothness = .5, phi = 1){
  # Simulates a mean-zero spatial field on the unit square
  #
  # Parameters:
  #  n - number of spatial locations
  #  range, smoothness, phi - parameters for Matern covariance function
  
  coords = matrix(runif(2*n), ncol=2)
  
  Sigma = Matern(d = as.matrix(dist(coords)), 
                 range = range, smoothness = smoothness, phi = phi)
  
  list(coords = coords,
       params = list(n=n, range=range, smoothness=smoothness, phi=phi),
       x = t(chol(Sigma)) %*%  rnorm(n))
}

# simulate data
x = simulate.field()

# configure gibbs sampler  
it = 100

# run sampler using default posteriors
post.samples = sFit(x = x$x, coords = x$coords, nSamples = it)

# build kriging grid
cseq = seq(0, 1, length.out = 10)
coords.krig = expand.grid(x = cseq, y = cseq)

# sample from posterior predictive distribution
burn = 75
samples.krig = sKrig(x$x, post.samples, coords.krig = coords.krig, burn = burn)

Named transformation functions

Description

Evaluates the named link function at the locations x.

Usage

tx(x, link, linkparams)

Arguments

Examples

bisque:::tx(0.5, 'logit', list(NA))

Derive parameters for building integration grids

Description

Note: w is defined on the transformed scale, but for convenience f is defined on the original scale.

Usage

wBuild(
  f,
  init,
  dim.theta2 = length(init),
  approx = "gaussian",
  link = rep("identity", length(init)),
  link.params = rep(list(NA), length(init)),
  optim.control = list(maxit = 5000, method = "BFGS"),
  ...
)

Arguments

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)))

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

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)))