maxp function

Maximum likelihood estimation

Find the maximum likelihood estimate for p, also equal probabilities

maxp(H, startp=NULL, give=FALSE, fcm=NULL, fcv=NULL, SMALL=1e-6, n=1,
   show=FALSE, justlikes=FALSE, ...)
maxplist(Hlist, startp=NULL, give=FALSE, fcm=NULL, fcv=NULL, SMALL=1e-6, ...)
maxp_single(H, startp=NULL, give=FALSE, fcm=NULL, fcv=NULL, SMALL=1e-6,
   maxtry=100, ...)
maxp_single2(H, startp=NULL, give=FALSE, fcm=NULL, fcv=NULL, SMALL=1e-6,
   maxtry=100, ...)
maxp_simplex(H, n=100, show=FALSE, give=FALSE, ...)
maxp_lsl(HLSL, startp = NULL, give = FALSE, fcm = NULL, fcv = NULL, SMALL=1e-6, ...)
equalp(H)

Arguments

• H: A hyper2 or hyper3 object
• Hlist: A list with elements all hyper2 objects
• HLSL: An lsl object
• startp: A vector of probabilities specifying the start-point for optimization; if a full unit-sum vector, then the fill-up value will be removed by indep() (except for maxp_lsl())
• give: Boolean, with default FALSE meaning to return just the evaluate (including fillup), and TRUE meaning to return the entire formal output of the optimization routine
• fcm,fcv: Further problem-specific constraints
• n: Number of start points to use
• show: Boolean, with TRUE meaning to show successive estimates
• justlikes: Boolean, with TRUE meaning to return just a vector of estimated likelihoods
• SMALL: Numerical minimum for probabilities
• maxtry: Integer specifying maximum number of times to try constrOptim() with slightly differing start points, to avoid a known bug which reports wmmin is not finite, bugzilla id 17703
• ...: Further arguments which maxp() passes to constrOptim()

Details

Function maxp() returns the maximum likelihood estimate for p, which has the unit sum constraint implemented.

Function maxplist() does the same but takes a list of hyper2 objects (for example, the output of ggrl()). Note that maxplist() does not have access to the gradient of the objective function, which makes it slow.

If function maxp() is given a suplist object it dispatches to maxplist().

Functions maxp_single() and maxp_single2() are helper functions which perform a single constrained optimization using base::constrOptim() or alabama::constrOptim.nl()

respectively. The functions should produce identical (or at least very similar) results. They are used by maxp() and maxp_simplex() which dispatch to either maxp_single() or maxp_single2() depending on the value of option use_alabama. If TRUE, they will use (experimental) maxp_single2(), otherwise (default) maxp_single(). Function maxp_single() is prone to the wmmin not finite bug [bugzilla id 17703] but on the other hand is a bit slower. I am not sure which one is better at this time.

Function maxp_simplex() is intended for complicated or flat likelihood functions where finding local maxima might be a problem. It repeatedly calls maxp_single(), starting from a different randomly chosen point in the simplex each time. This function does not take fcm or fcv arguments, it operates over the whole simplex (hence the name). Further arguments, ..., are passed to maxp_single().

The functions do not work for the masterchef_series6 likelihood function. These require a bespoke optimization as shown in the vignette.

Function equalp() returns the value of $p$ for which all elements are the same.

In functions maxp() etc, arguments fcm and fcv

implement linear constraints to be passed to constrOptim(). These constraints are in addition to the usual nonnegativity constraints and unit-sum constraint, and are added to the ui

and ci arguments of constrOptim() with rbind()

and c() respectively. The operative lines are in maxp_single():

UI <- rbind(diag(nrow = n - 1), -1, fcm)
    CI <- c(rep(SMALL, n - 1), -1 + SMALL, fcv)

where in UI, the first $n-1$ rows enforce nonnegativity of $p_i$, $1<=i<n$; row $n$ enforces nonnegativity of the fillup value $p_n$; and the remaining (optional) rows enforce additional linear constraints. Argument CI is a vector with corresponding elements.

Examples of their use are given in the icons vignette.

Author(s)

Robin K. S. Hankin

Note

In manpages elsewhere, n=2 is sometimes used. Previous advice was to use n=10 or greater in production work, but I now think this is overly cautious and n=1 is perfectly adequate unless the dimension of the problem is large.

The (bordered) Hessian is given by function hessian(), documented at gradient.Rd; use this to assess the sharpness of the maximum.

Function maxp() takes hyper2 or hyper3 objects but it does not currently work with lsl objects; use maxp_lsl().

The built-in datasets generally include a pre-calculated result of running maxp(); thus hyper2 object icons and icons_maxp are included in the same .rda file.

Function maxp() can trigger a known bug (bugzilla id 17703) which reports ‘wmmin is not finite’ . Setting option use_alabama to TRUE makes the package use a different optimization routine.

See Also

gradient,fillup

Examples

maxp(icons)

W <- hyper2(pnames=letters[1:5])
W1 <- ggrl(W, 'a', letters[2:3],'d')  # W1 is a suplist object
## Not run: maxp(W1)
  # takes a long time to maximize a suplist