sampling_nonlinear {multinomineq}

Posterior Sampling for Multinomial Models with Nonlinear Inequalities

Description

A Gibbs sampler that draws posterior samples of probability parameters conditional on a (possibly nonlinear) indicator function defining a restricted parameter space that is convex.

Usage

sampling_nonlinear(
  k,
  options,
  inside,
  prior = rep(1, sum(options)),
  M = 1000,
  start,
  burnin = 10,
  eps = 1e-06,
  progress = TRUE,
  cpu = 1
)

Arguments

k

vector of observed response frequencies.

options

number of observable categories/probabilities for each item type/multinomial distribution, e.g., c(3,2) for a ternary and binary item.

inside

an indicator function that takes a vector with probabilities p=c(p11,p12, p21,p22,...) (where the last probability for each multinomial is dropped) as input and returns 1 or TRUE if the order constraints are satisfied and 0 or FALSE otherwise (see details).

prior

a vector with two positive numbers defining the shape parameters of the beta prior distributions for each binomial rate parameter.

M

number of posterior samples drawn from the encompassing model

start

only relevant if steps is defined or cmin>0: a vector with starting values in the interior of the polytope. If missing, an approximate maximum-likelihood estimate is used.

burnin

number of burnin samples that are discarded. Can be chosen to be small if the maxmimum-a-posteriori estimate is used as the (default) starting value.

eps

precision of the bisection algorithm

progress

whether a progress bar should be shown (if cpu=1).

cpu

either the number of CPUs used for parallel sampling, or a parallel cluster (e.g., cl <- parallel::makeCluster(3)). All arguments of the function call are passed directly to each core, and thus the total number of samples is M*number_cpu.

Details

Inequality constraints are defined via an indicator function inside which returns inside(x)=1 (or 0) if the vector of free parameters x is inside (or outside) the model space. Since the vector x must include only free (!) parameters, the last probability for each multinomial must not be used in the function inside(x)!

Efficiency can be improved greatly if the indicator function is defined as C++ code via the function cppXPtr in the package RcppXPtrUtils (see below for examples). In this case, please keep in mind that indexing in C++ starts with 0,1,2... (not with 1,2,3,... as in R)!

For each parameter, the Gibbs sampler draws a sample from the conditional posterior distribution (a scaled, truncated beta). The conditional truncation boundaries are computed with a bisection algorithm. This requires that the restricted parameteter space defined by the indicator function is convex.

Examples

# two binomial success probabilities: x = c(x1, x2)
# restriction to a circle:
model <- function(x) {
  (x[1] - .50)^2 + (x[2] - .50)^2 <= .15
}

# draw prior samples
mcmc <- sampling_nonlinear(
  k = 0, options = c(2, 2),
  inside = model, M = 1000
)
head(mcmc)
plot(c(mcmc[, 1]), c(mcmc[, 2]), xlim = 0:1, ylim = 0:1)


##### Using a C++ indicator function (much faster)
cpp_code <- "SEXP inside(NumericVector x){
  return wrap( sum(pow(x-.50, 2)) <= .15);}"
# NOTE: Uses Rcpp sugar syntax (vectorized sum & pow)

# define function via C++ pointer:
model_cpp <- RcppXPtrUtils::cppXPtr(cpp_code)
mcmc <- sampling_nonlinear(
  k = 0, options = c(2, 2),
  inside = model_cpp
)
head(mcmc)
plot(c(mcmc[, 1]), c(mcmc[, 2]), xlim = 0:1, ylim = 0:1)