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.