Whereas the distribution function pnorMix is the trivial sum of
weighted normal probabilities (pnorm), its inverse,
qnorMix is computed numerically: For each p we search for
q such that pnorMix(obj, q) == p, i.e., \(f(q) = 0\)
for f(q) := pnorMix(obj, q) - p. This is a root finding
problem which can be solved by uniroot(f, lower,upper,*).
If length(p) <= 2 or method = "eachRoot", this happens
one for one for the sorted p's. Otherwise, we start by doing
this for the outermost non-trivial (\(0 < p < 1\)) values of p.
For method = "interQpspline" or "interpspline", we now compute
p. <- pnorMix(q., obj) for values q. which are a grid
of length l.interp in each interval \([q_j,q_{j+1}]\), where
\(q_j\) are the “X-extremes” plus (a sub sequence of length
n.mu.interp of) theordered mu[j]'s.
Then, we use montone inverse interpolation
(splinefun(q., p., method="monoH.FC")) plus
a few (maximally maxiter, typically one!) Newton steps.
The default, "interQpspline", additionally logit-transforms the
p. values to make the interpolation more linear.
This method is faster, particularly for large length(p).