count_binom {multinomineq}

Count How Many Samples Satisfy Linear Inequalities (Binomial)

Description

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.

Usage

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.

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

Value

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

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

See Also

bf_binom, count_multinom

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)