Get ML estimate for product-binomial/multinomial model with linear inequality constraints.
ml_binom(k, n, A, b, map, strategy, n.fit =3, start, progress =FALSE,...)ml_multinom(k, options, A, b, V, n.fit =3, start, progress =FALSE,...)
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.
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)
strategy: a list that defines the predictions of a strategy, seestrategy_multiattribute.
n.fit: number of calls to constrOptim .
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.
progress: whether a progress bar should be shown (if cpu=1).
...: further arguments passed to the function constrOptim. To ensure high accuracy, the number of maximum iterations should be sufficiently large (e.g., by setting control = list(maxit = 1e6, reltol=.Machine$double.eps^.6),outer.iterations = 1000.
options: number of observable categories/probabilities for each item type/multinomial distribution, e.g., c(3,2) for a ternary and binary item.
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.
Returns
the list returned by the optimizer constrOptim, including the input arguments (e.g., k, options, A, V, etc.).
If the Ab-representation was used, par provides the ML estimate for the probability vector θ.
If the V-representation was used, par provides the estimates for the (usually not identifiable) mixture weights α that define the convex hull of the vertices in V, while p provides the ML estimates for the probability parameters. Because the weights must sum to one, the α-parameter for the last row of the matrix V is dropped. If the unconstrained ML estimate is inside the convex hull, the mixture weights α are not estimated and replaced by missings (NA).
Details
First, it is checked whether the unconstrained maximum-likelihood estimator (e.g., for the binomial: k/n) is inside the constrained parameter space. Only if this is not the case, nonlinear optimization with convex linear-inequality constrained is used to estimate (A) the probability parameters θ
for the Ab-representation or (B) the mixture weights α for the V-representation.
Examples
# predicted linear order: p1 < p2 < p3 < .50# (cf. WADDprob in ?strategy_multiattribute)A <- matrix( c(1,-1,0,0,1,-1,0,0,1), ncol =3, byrow =TRUE)b <- c(0,0,.50)ml_binom(k = c(4,1,23), n =40, A, b)[1:2]ml_multinom( k = c(4,36,1,39,23,17), options = c(2,2,2), A, b
)[1:2]# probabilistic strategy: A,A,A,B [e1<e2<e3<e4<.50]strat <- list( pattern = c(-1,-2,-3,4), c =.5, ordered =TRUE, prior = c(1,1))ml_binom(c(7,3,1,19),20, strategy = strat)[1:2]# vertex representation (one prediction per row)V <- matrix(c(# strict weak orders0,1,0,1,0,1,# a < b < c1,0,0,1,0,1,# b < a < c0,1,0,1,1,0,# a < c < b0,1,1,0,1,0,# c < a < b1,0,1,0,1,0,# c < b < a1,0,1,0,0,1,# b < c < a0,0,0,1,0,1,# a ~ b < c0,1,0,0,1,0,# a ~ c < b1,0,1,0,0,0,# c ~ b < a0,1,0,1,0,0,# a < b ~ c1,0,0,0,0,1,# b < a ~ c0,0,1,0,1,0,# c < a ~ b0,0,0,0,0,0# a ~ b ~ c), byrow =TRUE, ncol =6)ml_multinom( k = c(4,1,5,1,9,0,7,2,1), n.fit =1, options = c(3,3,3), V = V
)