count_multinom: Count How Many Samples Satisfy Linear Inequalities (Multinomial)

Description

Draws prior/posterior samples for product-multinomial data and counts how many samples are inside the convex polytope defined by (1) the inequalities A*x <= b or (2) the convex hull over the vertices V.

Usage

count_multinom(
  k = 0,
  options,
  A,
  b,
  V,
  prior = rep(1, sum(options)),
  M = 5000,
  steps,
  start,
  cmin = 0,
  maxiter = 500,
  burnin = 5,
  progress = TRUE,
  cpu = 1
)

Value

a list with the elements

a matrix with the columns

count: number of samples in polytope / that satisfy order constraints
M: the total number of samples in each step
steps: the "steps" used to sample from the polytope (i.e., the row numbers of A that were considered stepwise)

with the attributes:

proportion: estimated probability that samples are in polytope
se: standard error of probability estimate

Arguments

k: the number of choices for each alternative ordered by item type (e.g. c(a1,a2,a3, b1,b2) for a ternary and a binary item type). The length of k must be equal to the sum of options. The default k=0 is equivalent to sampling from the prior.
options: number of observable categories/probabilities for each item type/multinomial distribution, e.g., c(3,2) for a ternary and binary item.
A: a matrix defining the convex polytope via A*x <= b. The columns of A do not include the last choice option per item type and thus the number of columns must be equal to sum(options-1) (e.g., the column order of A for k = c(a1,a2,a2, b1,b2) is c(a1,a2, b1)).
b: a vector of the same length as the number of rows of A.
V: a matrix of vertices (one per row) that define the polytope of admissible parameters as the convex hull over these points (if provided, A and b are ignored). Similar as for A, columns of V omit the last value for each multinomial condition (e.g., a1,a2,a3,b1,b2 becomes a1,a2,b1). Note that this method is comparatively slow since it solves linear-programming problems to test whether a point is inside a polytope (Fukuda, 2004) or to run the Gibbs sampler.
prior: the prior parameters of the Dirichlet-shape parameters. Must have the same length as k.
M: number of posterior samples drawn from the encompassing model
steps: an integer vector that indicates the row numbers at which the matrix A is split for a stepwise computation of the Bayes factor (see details). M can be a vector with the number of samples drawn in each step from the (partially) order-constrained models using Gibbs sampling. If cmin>0, samples are drawn for each step until count[i]>=cmin.
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.
cmin: if cmin>0: minimum number of counts per step in the automatic stepwise procedure. If steps is not defined, steps=c(1,2,3,4,...) by default.
maxiter: if cmin>0: maximum number of sampling runs in the automatic stepwise procedure.
burnin: number of burnin samples per step that are discarded. Since the maximum-likelihood estimate is used as a start value (which is updated for each step in the stepwise procedure in count_multinom), the burnin number can be smaller than in other MCMC applications.
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.

References

Hoijtink, H. (2011). Informative Hypotheses: Theory and Practice for Behavioral and Social Scientists. Boca Raton, FL: Chapman & Hall/CRC.

Examples

Run this code

### frequencies:
#           (a1,a2,a3, b1,b2,b3,b4)
k <- c(1, 5, 9, 5, 3, 7, 8)
options <- c(3, 4)

### linear order constraints
# a1 a2 < 1-a1-a2 <=> a1+2*a2 < 1)
# matrix A:
#             (a1,a2, b1,b2,b3)
A <- matrix(
  c(
    1, -1, 0, 0, 0,
    1, 2, 0, 0, 0,
    0, 0, 0, 1, -1,
    0, 0, 0, 0, 1
  ),
  ncol = sum(options - 1), byrow = TRUE
)
b <- c(0, 1, 0, .50)

# count prior and posterior samples and get BF
prior <- count_multinom(0, options, A, b, M = 2e4)
posterior <- count_multinom(k, options, A, b, M = 2e4)
count_to_bf(posterior, prior)

bf_multinom(k, options, A, b, M = 10000)
bf_multinom(k, options, A, b, cmin = 5000, M = 1000)