One-step-ahead (OSA) residuals via method = "cdf" are now available for many more distributions. During taping, pskewnorm(), pskewnorm2(), pskewt(), pskewt2() and pvm() integrate the density numerically with the AD-compatible integrate() of RTMB, which is therefore now required in version 2.0 or later. Added distribution functions for the beta-binomial (pbetabinom(), pztbetabinom(), pzibetabinom(), phbetabinom()) and the beta-negative binomial (pbnbinom(), pbnbinom2()), which sum the probability mass function, and for the Skellam distribution (pskellam()), which integrates it over the first mean. OSA residuals are supported for all of these and for the wrapped Cauchy. For the two circular distributions, the circle is cut at the fixed origin -pi, which keeps the residuals valid for hidden Markov and random effects models. Also fixed OSA residuals and simulation for dztbinom(), dztnbinom() and dztnbinom2(), and NaN OSA residuals for dcombinom() when an observation equals size.
Fixed two distribution functions that gave wrong values, and with them wrong OSA residuals: pzigamma() and pzigamma2() used the scale as a rate, so they were wrong whenever scale was not 1, and pzibinom() left out the binomial zeros at q = 0.
Fixed derivatives of several distribution functions, which matter for gradients, sdreport() and the Laplace approximation. pt.ad() had a zero derivative at 0, which gave wrong gradients in pt2(), pbct(), ptrunct() and ptrunct2() at q = mu, and made the derivative of dskewt() and pskewt() with respect to skew zero at skew = 0; the caution against starting skew at zero is therefore gone. ppowerexp(), ppowerexp2() and pbcpe() had NaN gradients at q = mu, pfoldnorm() could not be differentiated with respect to q, and pgenpois() not with respect to its parameters, so its OSA residuals failed with estimated parameters. The gamma-based distribution functions (pgamma2(), pzigamma(), pzigamma2(), pgengamma(), pinvgamma(), pinvchisq()) had non-finite second or third derivatives for small arguments, inherited from RTMB's pgamma(), which is now worked around.
dskewnorm(), dskewnorm2() and dexgauss() are now accurate far in the tails. They used log(1e-300 + pnorm(...)), which cut the log density off there and made its gradient zero, and now use pnorm(..., log.p = TRUE).
dinvgamma(), pinvgamma(), qinvgamma() and rinvgamma() can now be called with scale instead of rate, which previously failed. As in stats::dgamma(), giving both is an error.
Fixed the copula densities cgumbel() and cfrank(), which were wrong: the last factor of the Gumbel density and the sign in the denominator of the Frank density were incorrect, so neither integrated to one and likelihoods built with them through dcopula() were wrong. Both now agree with the copula package. The Frank density is also evaluated in a form that stays accurate for strong dependence. The copula distribution functions Cgumbel() and Cfrank(), used with ddcopula(), were correct and are unchanged.
Added two circular-linear copula constructors for use with dcopula(), joining a circular and a linear margin such as the turning angles and step lengths of an animal track. cjw() is the Johnson-Wehrly copula with any circular binding density, e.g. cjw(dvm, mu = 0, kappa = 2) or cjw(dwrpcauchy, mu = 0, rho = 0.5). Its dependence is a helix, so it is not symmetric in the sign of the angle. cfold() folds any of the linear copulas into a circular-linear copula that is symmetric, e.g. cfold(cgaussian(0.5)), and so links the straightness of a step to its length; this is the rectangular patchwork copula of Hodel and Fieberg (2022). Both allow for automatic differentiation.
Added the distribution function (pwrpcauchy()) and quantile function (qwrpcauchy()) of the wrapped Cauchy distribution, both in closed form. By default the circle is cut open at the antipode of the mean direction, mu - pi, the same origin as the default of pvm(), so pwrpcauchy(mu) is one half and angles outside (mu - pi, mu + pi] are wrapped onto that interval. A fixed origin can be set with from; this is needed whenever the distribution function is averaged over different values of mu, for example over hidden states or a random effect, because the average of distribution functions with different cuts is not a distribution function. pwrpcauchy() is differentiable, which makes the wrapped Cauchy usable as the turning-angle margin in copula models for step lengths and turning angles.
rwrpcauchy() now generates by inversion of pwrpcauchy() instead of calling circular::rwrappedcauchy(), so it recycles vector parameters. Random streams for a given seed differ from earlier versions, and with wrap = FALSE the angles now lie in [mu - pi, mu + pi] rather than [0, 2 * pi).
Added the half-t (dhalft()) and half-Cauchy (dhalfcauchy()) distributions, each with matching p, q and r functions. These are the standard weakly informative priors for the standard deviation of a hierarchical model, so they are aimed squarely at models fitted by the Laplace approximation. Both distribution functions are differentiable, and phalft() differentiates with respect to df as well as sigma, so the degrees of freedom can be estimated rather than fixed; one-step-ahead residuals via method = "cdf" are supported. The half-Cauchy is the half-t with df = 1, and the half-normal is dfoldnorm() with mu = 0, which the half-t approaches as df grows.
Added the Yule-Simon (dyules()) and Waring (dwaring()) distributions, two classical long-tailed count laws, each with matching p and r functions. Both are beta-geometric special cases of the beta-negative binomial, with size fixed at one, and that restriction is what gives them a closed-form distribution function where the general beta-negative binomial has none. One-step-ahead residuals via method = "cdf" are available. dyules() follows VGAM and is supported on the positive integers, while dwaring() follows the WARING family of gamlss.dist, starts at zero and has mu as its mean. The YULE family of gamlss.dist is dwaring(x, mu, mu).
Added the beta-negative binomial distribution (dbnbinom()) and its mean parameterisation (dbnbinom2()), each with a matching r function. It is to the negative binomial what the beta-binomial is to the binomial, and its extra beta layer gives a considerably heavier tail. In dbnbinom2(), which follows the BNB family of gamlss.dist, mu is exactly the mean; this is the more stable parameterisation to estimate in, because the original one has a long likelihood ridge along which size and shape2 trade off. The distribution function has no closed form and is summed from the probability mass function, see above.
Added the extreme value distributions: the generalised extreme value distribution (dgev()), the generalised Pareto distribution (dgpd()) and the Frechet distribution (dfrechet()), each with matching p, q and r functions. The first two cover their three shape regimes with a single expression rather than a branch on the sign of xi, so the derivative with respect to the shape is exact at xi = 0, which is the usual starting value when the shape is estimated. Densities and distribution functions are differentiable, so simulation and one-step-ahead residuals via method = "cdf" are supported. Unlike VGAM, evd and extraDistr, dgpd() returns 1 / sigma rather than zero at the threshold itself, matching stats::dexp() at zero.
RTMBdist no longer masks anything in stats. The AD-compatible replacements for stats::pt(), stats::plnorm(), stats::dgeom() and stats::pgeom() are exported as pt.ad(), plnorm.ad(), dgeom.ad() and pgeom.ad(), and are reached through internal S4 generics that dispatch on the argument classes: plain numeric input goes to the stats versions, AD variables to the .ad versions. Previously pt() was exported with a reduced argument list, so pt(q, df, lower.tail = FALSE) failed for anyone who had loaded the package. As a side effect lower.tail and log.p now also work for pt() under automatic differentiation.
Added the zero-inflated (dzigeom()), zero-truncated (dztgeom()) and hurdle (dhgeom()) geometric distributions. These build on RTMB's AD-compatible negative binomial with size = 1, so stats' own dgeom() and friends are left untouched.
Added the zero-inflated, zero-truncated and hurdle beta-binomial distributions (dzibetabinom(), dztbetabinom(), dhbetabinom()). Their distribution functions are summed from the probability mass function, see above.
The documentation of the continuous zero-inflated distributions now explains that, because the continuous part places no mass at zero, zeroprob is exactly the probability of a zero and zero-inflation coincides with a hurdle model; GAMLSS calls these zero-adjusted rather than zero-inflated.
Added the hurdle (zero-altered) count distributions: Poisson (dhpois()), binomial (dhbinom()), negative binomial (dhnbinom()) and its mean parameterisation (dhnbinom2()), each with matching p and r functions. In a hurdle distribution the probability of a zero is a free parameter and the positive counts follow the corresponding zero-truncated distribution. Unlike zero-inflation, which can only add zeros, a hurdle model allows zeroprob to be smaller than the Poisson would give on its own. The density and distribution function are both differentiable, so simulation and one-step-ahead residuals via method = "cdf" are supported.
Fixed pztnbinom() and pztnbinom2(), which returned NaN instead of 0 for quantiles below their support.
Added the Johnson SU distribution in both the original parameterisation (djsu(), pjsu(), qjsu(), rjsu()) and the moment parameterisation (djsu2(), pjsu2(), qjsu2(), rjsu2()), a four-parameter distribution on the real line covering a wide range of skewness and kurtosis. In djsu2() the location and scale arguments are exactly the mean and standard deviation. The density and distribution function are both differentiable, so simulation and one-step-ahead residuals are supported. Unlike gamlss.dist, the reparameterisation stays finite for very large tau, where the distribution approaches the normal.
Added references to the primary source for each distribution derived from gamlss.dist, and cross-links between related families. pgenpois() and friends previously had no references at all.
Removed the dependency on gamlss.dist, which is scheduled for archival on CRAN. The quantile and random generation functions of the Box-Cox Cole-Green (qbccg(), rbccg()), Box-Cox t (qbct(), rbct()), Box-Cox power exponential (qbcpe(), rbcpe()), power exponential (qpowerexp(), rpowerexp(), qpowerexp2(), rpowerexp2()), Pareto (qpareto(), rpareto()), generalised Poisson (pgenpois(), qgenpois(), rgenpois()) and exponentially modified Gaussian (qexgauss()) distributions are now implemented natively. Results are unchanged except for the fixes below.
Fixed argument recycling in qbccg(), qbct(), qbcpe(), qgenpois() and pgenpois(). Evaluating a single quantile against vectorised parameters previously collapsed the result to length one, and a parameter vector shorter than x silently truncated it. These functions now return one value per recycled argument tuple.
Fixed qbct() and rbct(), which failed with an error when mu, sigma or tau was supplied as a vector.
qpareto() now honours lower.tail, which was previously accepted but ignored.
log.p = TRUE now works in qpareto(), qbct() and qgenpois(). These previously validated p before transforming it back from the log scale, so the argument could not be used.
qexgauss() is substantially more accurate. It inverts pexgauss() with a tighter convergence tolerance, reducing the round-trip error |p(q(p)) - p| from roughly 1e-6 to roughly 1e-14.
pgenpois() no longer falls back to the Poisson distribution for phi < 1e-4 and evaluates the generalised Poisson distribution function across the whole parameter range.
pbetaprime(), pinvchisq() and pinvgamma() are now differentiable, so one-step-ahead residuals via method = "cdf" are available for the beta prime, inverse chi-squared and inverse gamma distributions.
Added the Bell distribution (dbell(), pbell(), qbell(), rbell()) and its mean parameterisation (dbell2(), pbell2(), qbell2(), rbell2()), a one-parameter distribution for overdispersed counts. Log Bell numbers are used throughout and cached, so the density stays finite well past the point at which the Bell numbers themselves overflow double precision. One-step-ahead residuals are supported via method = "cdf".
Added lambertW(), an AD-compatible implementation of the principal branch of the Lambert W function. Derivatives of all orders are exact, so it can be used in models fitted by Laplace approximation.
Added the Conway-Maxwell-binomial distribution (dcombinom(), pcombinom(), qcombinom(), rcombinom()), a binomial generalisation with an additional dispersion parameter nu. The density and distribution function are differentiable with respect to x, so one-step-ahead residuals are supported.