count_binom() R function from [multinomineq]

Count How Many Samples Satisfy Linear Inequalities (Binomial)

Draws prior/posterior samples for product-binomial 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.


count_binom(
  k,
  n,
  A,
  b,
  V,
  map,
  prior = c(1, 1),
  M = 10000,
  steps,
  start,
  cmin = 0,
  maxiter = 500,
  burnin = 5,
  progress = TRUE,
  cpu = 1
)

Arguments

k: vector of observed response frequencies.
n: the number of choices per item type. If k=n=0, Bayesian inference is relies on the prior distribution only.
A: a matrix with one row for each linear inequality constraint and one column for each of the free parameters. The parameter space is defined as all probabilities x that fulfill the order constraints A*x <= b.
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.
map: optional: numeric vector of the same length as k with integers mapping the frequencies k to the free parameters/columns of A/V, thereby allowing for equality constraints (e.g., map=c(1,1,2,2)). Reversed probabilities 1-p are coded by negative integers. Guessing probabilities of .50 are encoded by zeros. The default assumes different parameters for each item type: map=1:ncol(A)
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
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.

Returns

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
const_map: logarithm of the binomial constants that have to be considered due to equality constraints

Details

Counts the number of random samples drawn from beta distributions that satisfy the convex, linear-inequalitiy constraints. The function is useful to compute the encompassing Bayes factor for testing inequality-constrained models (see bf_binom; Hojtink, 2011).

The stepwise computation of the Bayes factor proceeds as follows: If the steps are defined as steps=c(5,10), the BF is computed in three steps by comparing: (1) Unconstrained model vs. inequalities in A[1:5,]; (2) use posterior based on inequalities in A[1:5,] and check inequalities A[6:10,]; (3) sample from A[1:10,] and check inequalities in A[11:nrow(A),] (i.e., all inequalities).

Examples


### a set of linear order constraints:
### x1 < x2 < .... < x6 < .50
A <- matrix(
  c(
    1, -1, 0, 0, 0, 0,
    0, 1, -1, 0, 0, 0,
    0, 0, 1, -1, 0, 0,
    0, 0, 0, 1, -1, 0,
    0, 0, 0, 0, 1, -1,
    0, 0, 0, 0, 0, 1
  ),
  ncol = 6, byrow = TRUE
)
b <- c(0, 0, 0, 0, 0, .50)

### check whether a specific vector is inside the polytope:
A %*% c(.05, .1, .12, .16, .19, .23) <= b

### observed frequencies and number of observations:
k <- c(0, 3, 2, 5, 3, 7)
n <- rep(10, 6)

### count prior samples and compare to analytical result
prior <- count_binom(0, 0, A, b, M = 5000, steps = 1:4)
prior # to get the proportion: attr(prior, "proportion")
(.50)^6 / factorial(6)

### count posterior samples + get Bayes factor
posterior <- count_binom(k, n, A, b, M = 2000, steps = 1:4)
count_to_bf(posterior, prior)

### automatic stepwise algorithm
prior <- count_binom(0, 0, A, b, M = 500, cmin = 200)
posterior <- count_binom(k, n, A, b, M = 500, cmin = 200)
count_to_bf(posterior, prior)

References

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

Fukuda, K. (2004). Is there an efficient way of determining whether a given point q is in the convex hull of a given finite set S of points in Rd? Retrieved from https://www.cs.mcgill.ca/~fukuda/soft/polyfaq/node22.html

count_binom function