Added the Poisson-inverse Gaussian (dpig()) and Sichel (dsichel()) distributions, each with matching p and r functions. The Poisson-inverse Gaussian is the standard alternative to the negative binomial for overdispersed counts: a Poisson whose mean follows an inverse Gaussian, parameterised as the PIG family of gamlss.dist, so that mu is the mean and the variance is mu + sigma * mu^2, the same as for dnbinom2() with size = 1 / sigma. Compared with the negative binomial it has a heavier right tail and more zeros. The Sichel distribution, the SICHEL family of gamlss.dist, mixes the Poisson over a generalised inverse Gaussian instead and adds a shape parameter nu; the Poisson-inverse Gaussian is its special case nu = -1/2. Both mass functions involve Bessel functions whose order grows with the count, which besselK() cannot evaluate for large counts and strong overdispersion because it overflows; for orders above 50, the uniform asymptotic expansion of the Bessel function is used instead, so the mass functions and their derivatives stay finite and accurate for any count. nu can be estimated, as the derivative with respect to the order is available. The distribution functions are summed from the mass functions, so one-step-ahead residuals via method = "cdf" are supported.
Added the wrapped normal (dwrpnorm()) and cardioid (dcard()) distributions, each with matching p and r functions, joining the von Mises and wrapped Cauchy as turning-angle distributions for movement models and as binding densities for cjw(). The wrapped normal density is a sum over the windings of the normal density for sigma < 2 and its Fourier series otherwise, each truncated where the remaining terms are below 1e-20. The distribution functions are in closed form and cut the circle open where pwrpcauchy() does, at mu - pi or at from, and one-step-ahead residuals via method = "cdf" cut it at -pi, as for the other circular distributions.
Added the generalised beta distribution of the second kind (dgb2()) and its special cases, the Dagum (ddagum()), Singh-Maddala or Burr XII (dsinmad()) and Lomax or Pareto II (dlomax()) distributions, each with matching p, q and r functions. These are the standard models for income, wealth and claim sizes. The distribution functions are incomplete beta functions, computed with the upper tail directly rather than as 1 - p, and all of them are differentiable, so one-step-ahead residuals via method = "cdf" are supported. The log-logistic distribution dllogis() is the GB2 with shape2 = shape3 = 1.
Added the Birnbaum-Saunders or fatigue life distribution (dbisa()), with matching p, q and r functions, the standard model for time to failure under cumulative damage. Density, distribution and quantile function are in closed form and differentiable, and both tails of the distribution function are accurate; the quantile function is written as scale * exp(2 * asinh(shape * z / 2)), which avoids the cancellation of the textbook formula for small probabilities.
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 wrong gradients and Hessians of the distribution functions built on the incomplete beta function: pbeta2(), qbeta2(), pbetaprime(), qbetaprime(), pzibeta(), pzoibeta(), poibeta(), pnbinom2(), pztnbinom(), the new GB2 family and the families built on these. RTMB's pbeta() and qbeta() (version 2.0) have wrong, often zero, derivatives in the shape parameters when a shape is an integer or half-integer, as for integer counts in pnbinom2() or integer starting values, and inaccurate second derivatives close to such values. This is fixed in the development version of TMB; until it reaches CRAN, the shapes are moved slightly off these values and the results extrapolated back, which keeps values accurate to about 1e-11.
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.
rskewt() and rskewt2() now recycle vector parameters. They called sn::rst(), which is not vectorised over its parameters, and failed for example when simulating from a hidden Markov model with state-dependent parameters. The same held for pskewt(), qskewt(), pskewt2() and qskewt2() with vector skew or df, which now call sn::pst() and sn::qst() once for each distinct pair of skew and df. rskewt() and rskewnorm() (and with them rskewt2() and rskewnorm2()) now draw from the stochastic representation of the skew normal and skew t distributions instead of calling the sn package, which is also much faster for vector parameters. Random streams for a given seed differ from earlier versions.
pskewt(), pskewt2(), pskewnorm() and pskewnorm2() are now accurate in both tails. Outside of AD they used the sn package throughout, whose values are inaccurate in the tails: the upper tail was computed as 1 - p, which loses all precision there, sn::psn() is far off in the light tail (e.g. 7e-22 instead of 9e-31 at z = -5, alpha = 2), and sn::pst() is off by over 1% beyond 30 scale units for non-integer df. Where the requested tail probability is below 1e-5 or q is that far out, the tail is now integrated numerically, as is already done during taping, which also makes log.p = TRUE accurate where the probability underflows.
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 three examples to the Worked Examples vignette: the five largest annual sea levels at Venice with a latent random-walk trend (dgev(), pgev()), the ETAS model for earthquake aftershocks (dlomax(), plomax(), dgpd()), and pooling opinion polls with house effects (dbetabinom()).
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.