Mainly a bug-fix release, written for users of 0.1.1; the default method of
rmhn() also chooses its sampler by re-derived thresholds (see Seeded draws
below and Performance).
0.1.1 returns wrong values without a warning in several places, some of them
at ordinary parameter values: the skewness and kurtosis at a large shape
(mhn_kurtosis(1000, 1, 0) is -0.627 where the true value is 7.5e-7), the
variance for a shape close to 1, the distribution function under a strong
negative tilt when the mode or the quantile lies below about 1.5e-8, the
logarithm of an upper tail close to one and the quantile of such a tail,
and the density, distribution function and quantile for a shape or a tilt
within 1.5e-8 of alpha = 1 or gamma = 0 but not equal to it. A
few of these came with 0.1.1 itself, and 0.1.0 is right there. The true values
quoted below come from independent references: arbitrary-precision
computations, or the Gamma and normal laws the distribution approaches in the
limit. Each fix carries a regression test at the values that failed.
Deterministic functions. On the grid of 330 parameter triples described
under 0.1.1 below, in the 130 cells with |gamma| / sqrt(beta) at most 10,
dmhn(), pmhn(), mhn_mean() and mhn_var() agree with 0.1.1 to twelve
significant digits in 128 cells and differ by at most 5.3e-11 relative in the
other two; qmhn() at p = 0.001, 0.1, 0.5, 0.9 and 0.999 differs in 13
cells, by at most 5.9e-9; mhn_mode() does not change. mhn_skewness() and
mhn_kurtosis() change by up to 2.2e-5 and 2.8e-4 relative, for
alpha <= 1/2 at gamma / sqrt(beta) = 10, where 0.1.1 is the inaccurate one
(see below). Outside these cells values change by more, and those changes are
the fixes listed below.
Seeded draws. For a fixed seed rmhn() returns the draws of 0.1.1 over
most of the parameter space wherever it uses the sampler 0.1.1 used. The
default method now uses RTDR in parts of the space where 0.1.1 used Sun et
al.'s Algorithm 1 or 3 (see Performance), and there its draws follow a
different sequence. On the 330-cell grid (200 draws per cell under
set.seed(20260903), default method), 300 cells give the same draws and 30
a different sequence: 3 (alpha < 1/2 at a tilt of 1e7) where RTDR itself
changed, and 27 where the default method changed sampler (12 with
gamma > 0 and alpha > 1 at a tilt up to 2 sqrt(alpha), 15 with
gamma < 0 and alpha >= 10). In more detail:
gamma / sqrt(beta) of about 5e6
for alpha < 1/2, 1.6e7 for 1/2 <= alpha < 1 and 7.7e6 for alpha > 1
(6.1e6 at a shape of 1e12), and at every negative tilt checked (down
to -1e8), both on the default path and under method = "rtdr", at shapes
up to 1e12 (apart from a rare change in the trailing digits at an extreme
scale); the
square-root-Gamma case gamma = 0 for shapes from 0.05 to 1e16; and the
truncated normal of alpha = 1 up to a tilt of about 1e8, or at any tilt
when beta = 1.gamma > 0, alpha > 1), in the last bit or two, for shapes
below 1e6 at every beta checked (1e-20 to 1e20), as does the truncated
normal of alpha = 1 beyond a tilt of 1e8 with beta != 1; and
Algorithm 3, by at most about 1e-12 relative for beta <= 1, a few 1e-12
for beta up to 10, 1e-9 at beta = 1e6 and 1e-7 at beta = 1e10, and
at beta = 1e20, where the search of 0.1.1 for the matching point had
lost precision, by up to 4e-3 relative (0.005 standard deviations).
Algorithm 1 supplies every draw with gamma > 0 and alpha > 1 under
method = "sun", and those of the default method there at few draws per
setup or a tilt above 2 sqrt(alpha); Algorithm 3 supplies every draw with
gamma < 0 under method = "sun", and those of the default method for
gamma < 0 with alpha < 1 or alpha >= 1000 at few draws per setup,
past a tilt of 1e150 and where RTDR cannot place the law (?rmhn gives
the rule). An
accept-reject decision that fell within such a change would shift the
sequence from there on; none did in the draws compared.gamma / sqrt(beta)
of about 52.5 to 55, from the first draw whose log-axis value would pass
the largest double and that the corrected test then rejects; such a draw
has a chance of about exp(-1.8e308 alpha) per proposal from the left
tail (e^-4 at 2^-1022, 1.6e-8 at 1e-307). Where every draw is 0 the values
stay and the stream after them can move (see the RTDR entry under the
random-generation fixes).1 / (1 + w), decided by one uniform variate, and
otherwise drawn as from the shape 2^-60 (see the subnormal-shape entry
below); below a tilt of about 52 (53 at the smallest shapes) every draw
is 0 and uses no random number, so that the draws after it in the stream
change too.alpha > 1, at the large
shapes where it builds its envelope relative to the mode (from about
2.4e12 to 3.1e12 at moderate tilts, 2.3e13 near a tilt of -alpha / e and
4.5e13 at any tilt); gamma = 0, and Algorithm 3, at a shape below
about 0.05, from the first Gamma variate that underflows to zero onwards
(at gamma = 0 this happens with probability about 1.5e-5 per draw at
alpha = 0.03 and 1e-8 at alpha = 0.05); Algorithm 1 at a shape of
about 1e6 or more with beta != 1; shapes above about 1e14 under a
negative tilt, above about 1e16 under a positive one and above 2^60
(about 1.2e18) at gamma = 0; and every parameter set at which 0.1.1
returned NaN, stopped with an error, or sampled the wrong law (see
below).gamma < 0 with
1 < alpha < 1000 at fewer than 25 draws per setup or with
alpha >= 10, and with alpha >= 1000 from 10, 25 or 50 draws per
setup, and at any number where RTDR's setup starts next to the mode;
gamma < 0 with 0.1 <= alpha < 1 at 5 to 24 draws per setup and with
alpha < 0.1 at 25 to 99; and gamma > 0 with alpha > 1 at a tilt up
to 2 sqrt(alpha) from 5 to 1000 draws per setup, by shape and tilt.
This includes the single draws of a Gibbs sampler with gamma < 0 and
1 < alpha < 1000. With recycled parameter vectors the draws per setup
are also counted differently (see Performance), which can change the
sampler, and with it the sequence.Invalid arguments. dmhn(), pmhn(), qmhn() and rmhn() now treat a
parameter outside the parameter space, NA and a zero-length argument as the
base R distribution functions do, element by element, where 0.1.1 stopped
the whole call; and a factor or a string is an error, where 0.1.1 used it
as a number (see Interface). Valid arguments give the values and the draws
they gave before.
mhn_skewness() and mhn_kurtosis() were wrong at a large shape under a
non-positive tilt, because their quadrature range ended where the kernel had
fallen by only a few units of log density: mhn_kurtosis(1000, 1, 0)
returned -0.627 and mhn_skewness(1000, 1, 0) 0.00106, where the true
values are 7.5e-7 and 0.0224; at alpha = 1e4 the kurtosis was -1.14, and
mhn_kurtosis(200, 1, -2) had the wrong sign (-0.0040 against 0.00044).
For alpha < 1 with a positive tilt strong enough to give the density an
interior peak, the quadrature panel next to the origin contained the
singularity of x^(alpha - 1): mhn_skewness(0.01, 1, 3) and
mhn_kurtosis(0.01, 1, 3) returned 1.20 and 0.46 against 3.47 and 12.15,
and mhn_kurtosis(0.3, 1, sqrt(10)) was -0.55609 against -0.55730.
0.1.0 was right at these values.
The normalising constant was taken from its alpha = 1 or alpha = 2
closed form for any shape within sqrt(.Machine$double.eps) of those
values, which put it out of step with the rest of the computation:
mhn_var(1 + 1e-9, 1, -20) was 7.4e-7 relative off,
mhn_var(1 + 1.4e-8, 1, -1e6) returned 1.0e-7 where the variance is
1.0e-12, and the density at alpha = 2 + 1e-8, gamma = 5 integrated to
1 + 9.6e-9. The closed forms are now used only at alpha = 1 and 2
exactly. mhn_mode() likewise treated a shape near 1 as 1:
mhn_mode(1 + 1e-9, 1, -1) was 0, where the mode is 1.0e-9.
mhn_var() used the closed form of Sun et al. (2023, Lemma 2c) where its
cancellation was larger than its guard allowed for, and returned
variances above the bound 1/(2 beta) that Sun et al. (2023, Lemma 4c)
prove for alpha >= 1: 0.50127 at (2, 1, 900) against a true 0.4999988,
6.17 at (1e8, 1, 1) and 1074 at (1e9, 1, -1), where the variance is
about 0.25.
Each moment helper now bounds the error of its closed form and integrates
directly where that bound exceeds the working tolerance.
mhn_mean(alpha, beta, 0) lost accuracy at a large shape: at alpha = 1e12
it was 1.9e-4 low, and at alpha = 1e15 21.7 times too large. It is now
accurate to the rounding of its logarithm, about log(alpha) / 2: within
4e-15 relative up to alpha = 1e15 and 4e-14 at any shape up to the
largest double.
mhn_mode() returned the inflection point gamma / (4 beta) at the
boundary alpha = 1 - gamma^2 / (8 beta) of Sun et al. (2023, Lemma 3c),
where the density has no interior maximum: mhn_mode(0.5, 1, 2) was 0.5.
It now returns NA there, as for any alpha < 1 without an interior mode.
(Lemma 3c states its condition with <= and places a local maximum there;
at equality the stationary points merge into an inflection point, so the
package departs from the lemma at this boundary.)
mhn_kurtosis(2, 1e300, 1) was NaN; it is now 0.2451.
For alpha just below 1 and a small positive tilt, just past the tilt
gamma / sqrt(beta) = sqrt(8 (1 - alpha)) at which the density acquires
its interior peak, mhn_skewness() and mhn_kurtosis() were up to 1e-7
off: mhn_kurtosis(0.99998, 1, 0.013) returned 0.85932105 where the true
value is 0.85932098. The quadrature now takes the bend of x^(alpha - 1)
right of the peak on its own, and on a grid of shapes from 1e-6 to
1 - 1e-8 with an interior peak they are within 1.2e-11 of
arbitrary-precision references.
At alpha = 1 under a positive tilt gamma / sqrt(beta) from about 11.1
on, mhn_mean() returned values just below the mode gamma / (2 beta) and
mhn_var() values just above 1/(2 beta), against E(X) >= mode, which
the truncated normal satisfies and Sun et al. (2023, Lemma 4a) prove for
alpha > 1, and the bound Var(X) <= 1/(2 beta) of their Lemma 4c: up to a
tilt of 12 by up to 6e-15 of the mean and 2.2e-13, and beyond by more (8e-12
of the mean at tilts up to 1000; for the variance see the item on the
Lemma 2c form above). 0.1.1 did the same for a shape within 1.5e-8 of 1,
which it treated as 1. Where the closed forms of Lemmas 2a and 2c would
still cross these bounds next to alpha = 1, the quadrature is now taken
instead, and for alpha >= 1 a variance from the quadrature is capped at
1/(2 beta), which is then the nearer value.
Under a strong negative tilt mhn_mean() could exceed the upper bound of
Lemma 4a. At a shape of 6 or more it was off by up to 2.8e-9
(mhn_mean(3e6, 1, -35) returned 1216.02602953 where the mean is
1216.02602617), and mhn_mean(1e5, 1, -1e9) returned 1.0000000000387e-4,
where the mean and the bound both round to 0.9999999999998e-4. There the
mean now comes from the quadrature, which keeps to the bound to within an
ulp and is within 3e-15 of arbitrary-precision references for shapes from
6 to 3e6. Below a shape of 6 it was off by up to 3.6e-12 for tilts down to
-1e14, and above the bound at about a third of the points checked. It is
still the Lemma 2a ratio there, now within 5e-14 of the references down to
that tilt, and where that error would carry it above the bound, from a
tilt of about -1e7 on, the bound is returned instead, which is closer to
the mean.
For a small shape under a positive tilt the normalising constant came
from Sun et al.'s Supplementary Lemma 10 series, whose second term was
formed from (alpha / 2 + 1) - 1: that keeps only the leading bits of
alpha / 2, and none below 1.1e-16, and every later term of the even half
of the series inherited the error, about 2.2e-16 / alpha relative.
Wherever the peak carries the mass, from a tilt of about 12 at such a
shape, the density and the moments were off with it: at alpha = 1e-16
and a tilt from 15 up the density was twice its value and
mhn_mean(1e-16, 1, 20) 19.90 against 9.949; at 1e-15 and 1e-14 the
density was 6 and 1 percent off at a tilt of 20, and at 1e-12 the mean
4.4e-5 off. The term is now formed as alpha / 2 + (k - 1), and past a
tilt of 20 the constant comes from quadrature (see below).
dmhn(), pmhn() and qmhn() treated a shape within
sqrt(.Machine$double.eps) (1.5e-8) of 1 as 1 and a tilt
|gamma| / sqrt(beta) below it as 0, and answered for the truncated normal
or the square root of a Gamma law; the normalising constant, and with it
the moments, did the same for the tilt. The law requested differs from
those by about (alpha - 1) (log(sqrt(beta) x) - E log U) or
gamma / sqrt(beta) (sqrt(beta) x - E U) in the log density,
U = sqrt(beta) X: dmhn(1e-300, 1 + 1.4e-8, 1, -1) was 9.7e-6 relative
off, pmhn(6, 2.5, 1, 1.4e-8, lower.tail = FALSE) 7.1e-8 and
mhn_mean(2.5, 1, 1.4e-8) 3.1e-9. The closed forms are now used at
alpha = 1 and gamma = 0 exactly, and the general forms, which keep
their accuracy up to those values, everywhere else. On 126 triples next to
them the density, both tails and the quantile agree with an
arbitrary-precision reference to within 3.3e-11 up to a shape of 50
(0.1.1: up to 9.7e-6 in the density), and above it to within the ordinary
accuracy of the general forms, 9e-11 at alpha = 1e5, which there exceeds
the error of taking a tilt below about 1e-12 as 0. rmhn() still draws
from those laws within the same window, as its help page says: the law
drawn is within a total variation
distance of 8.6e-9 of the one requested, and the draws for a fixed seed do
not change. (The general samplers draw the exact law there, but of Sun et
al.'s samplers only Algorithm 2, which the package does not implement,
covers alpha < 1 with gamma > 0, so method = "sun" would stop next to
alpha = 1 from below under a positive tilt, and next to gamma = 0 from
above for any alpha < 1.)
For alpha above about 1.0e7 with a nonzero tilt, the normalising constant
lost a factor of 2 and dmhn() returned twice the density:
dmhn(mhn_mode(1.5e7, 1, 1), 1.5e7, 1, 1) was 1.5958, where the density is
0.79785. pmhn() and qmhn() were wrong in the same range:
pmhn(mhn_mode(1.5e7, 1, 1), 1.5e7, 1, 1) was 1 and
pmhn(mhn_mode(1e8, 1, 1), 1e8, 1, 1) 0, where both are 0.5, and
qmhn(0.5, 1e8, 1, 1) was 7071.66, where the median is 7071.32.
Under a negative tilt at a large shape the quadrature for the normalising
constant could miss part of the mass: at (2795.013, 1, -9.112) the density
was 9.6e-5 relative too low. For a shape below 1 its logarithm was off by
up to 4.4e-10, at (0.7148, 1, -0.6437), and by more than 1e-12 at one in
a hundred of 9.4e4 such points checked; the density and the distribution
function carry that error. Under a negative tilt it is now within
4.4 eps (1 + |log Psi|) of an arbitrary-precision reference at 1.27e5
points with shapes from 1e-4 to 2000 and gamma / sqrt(beta) from -1e5
to -1e-4.
For gamma / sqrt(beta) above 20, and for alpha above 100, the
normalising constant now comes from quadrature, which there is both faster
and more accurate than the Lemma 10 series: dmhn(500, 2.5, 1, 1000) was
4.8e-9 relative off and took more than 10 ms; it now takes a few
microseconds.
The quadratures held the error of each panel to a tolerance that grew with
the panel's width, so their accuracy and cost depended on beta, on which
the shape of the law does not depend. At a large beta the tolerance could
not be met: mhn_skewness(1.5, 1e12, -3e6) took 38 ms, against about 45
microseconds for the same law at (1.5, 1, -3). Both now take about 9
microseconds.
A subnormal shape, below the smallest normal double 2.2e-308, was not
served correctly. There half of alpha is not always a double (half of
3 * 2^-1074 rounds to 2 * 2^-1074, and half of 2^-1074 to 0), and
1 / alpha passes the largest double below 5.6e-309. The series and the
square-root-Gamma case that start from alpha / 2, and the quadratures
that hold the mass by the origin, about 1 / alpha, went wrong with
them: dmhn(1, 2^-1074, 1, 1, log = TRUE) returned Inf and
pmhn(0.5, 2^-1074, 1, 1) NA; at 3 * 2^-1074 the log density was
-743.054 against -743.341 (a third too large as a density),
pmhn(1e-300, 3 * 2^-1074, 1, 1) was 0 against 1, mhn_mean(3 * 2^-1074, 1, 12) 1.510e-307 against 1.133e-307 and
qmhn(1e-320, 3 * 2^-1074, 1, 1, lower.tail = FALSE) 4.9e-324 against
2.56e-293; and mhn_skewness(2 * 2^-1074, 1, -30) was -Inf against
6.34e161. At such a shape t^alpha is 1 to double precision at every
positive double t, and the law takes a form in which alpha only
multiplies: the density is
alpha exp(-beta x^2 + gamma x) / (x (1 + w)), the upper tail
alpha J(sqrt(beta) x) / (1 + w) with
J(u) = int_u^Inf exp(-t^2 + z t) / t dt, the lower tail
(1 + alpha (K0(z) - J(sqrt(beta) x))) / (1 + w), and the moments
E X^k = alpha A_k / ((1 + w) beta^(k/2)) with
A_k = int_0^Inf t^(k-1) exp(-t^2 + z t) dt, z = gamma / sqrt(beta),
w = alpha K0(z) and
K0(z) = int_0^1 (exp(-t^2 + z t) - 1) / t dt + int_1^Inf exp(-t^2 + z t) / t dt.
1 / (1 + w) is the probability that X lies below the smallest positive
double; it is 1 in double precision below a tilt of about 52 (53 at the
smallest shapes), and for a
large tilt w is the mass at the peak of the kernel against the mass by
the origin. dmhn(), pmhn(), qmhn(), rmhn() and the four moment
functions now use these forms for a subnormal shape up to a tilt between
about 55 and 56.3, as alpha falls. Above it the mass by the origin is at
most 2^-64 (5.4e-20) of the whole, the density, the upper tail and the
moments are those of the ordinary paths, and the draws follow the same
law; the lower tail far below the peak, of which that mass is more than
2^-64, and its quantiles still come from these forms up to one half
(pmhn(1, 2^-1074, 1, 57, log.p = TRUE) was 0 against -65.03). rmhn()
draws 0 with probability 1 / (1 + w): at 2^-1074 and a tilt of
54.669, where that is 1/2, it drew 0 for 73 percent of the draws. While the mass
at the peak is small (w well below 1, up to a tilt of about 53 to 54.5
as alpha falls), the skewness and the excess kurtosis grow like
1 / sqrt(alpha) and 1 / alpha, by factors that depend on the tilt,
and the kurtosis is Inf where it passes the largest double: below a
shape of about 1.1e-308 at z = 0, at every subnormal shape for
z = -30, and at none for z = 12. mhn_mode() was already right
there.
At gamma = 0, qmhn() with log.p = TRUE lost the far tail opposite
to the one given. That case hands the probability to R's qgamma()
(X^2 follows a Gamma law), which returns a wrong quantile for a log
probability close to 0 at some shapes, of either tail:
qmhn(log1p(-1e-300), 1e-10, 1, 0, log.p = TRUE) returned 10.92, where
the quantile is 25.70, and at alpha = 10 the quantile was 0.18 to 0.33
percent off; qmhn(-1e-300, 500, 1, 0, lower.tail = FALSE, log.p = TRUE)
returned 2.4267, where the quantile is 2.4561, and at alpha from about
200 to 3e4 it was up to 8.4 percent off. A log probability above
log(1/2) is now passed to qgamma() as the logarithm of the other tail.
The lower-tail quadrature of pmhn(), which serves a strong negative tilt,
normalised its kernel at a point floored at sqrt(.Machine$double.eps), so
it was wrong whenever the mode or q lay below about 1.5e-8:
pmhn(2e-9, 3, 1, -1e8) returned 0.0175 and pmhn(1e-9, 3, 1, -5e6)
4.3e-6, where the law is Gamma(3, rate -gamma) to within 1e-11 and the
true values are 0.00115 and 2.1e-8. qmhn() inherited the error:
qmhn(0.2, 3, 1, -3e8) returned 8.56e-9 against 5.12e-9, and
qmhn(0.5, 3, 1, -3e8) failed with a warning.
Far from the peak the lower-tail quadrature could fall back towards
P(X <= mode) and stop being monotone: pmhn(c(0.01, 1, 6), 2.5, 1, -3e4)
was 1 1 0.300, and pmhn(2.47e8, 0.7, 1, -20) was 0.9595 where it is 1.
pmhn(q, lower.tail = FALSE, log.p = TRUE) took the logarithm of the
upper-tail quadrature, so where P(X > q) is close to one its logarithm,
about -F(q), kept only the absolute accuracy of the quadrature: at the
1e-10 quantile of (2, 1, 1) it returned -9.9999953e-11 for -1e-10, and at
that of (3, 1, 10) -1.000195e-10. That absolute error is about eps
times the size of the terms the logarithm is assembled from,
|log Psi| + (alpha/2) |log beta| + m (1 + |log v*|), with
gamma^2 / (4 beta) removed from log Psi for gamma / sqrt(beta) above
1000 (see ?pmhn). Wherever it could pass 1e-9 of F, below the median,
the logarithm is now formed as log1p(-F) from the lower tail: at
parameters of order one below F of about 1e-5, and at every such q
where the terms are largest (at (2e4, 1e300, 5e151) the logarithm of the
quadrature is 2.6e-8 off at F = 0.05). Elsewhere it stays the logarithm
of the quadrature, whose relative error there is held below 1e-9 by the
error bound it is gated on: four times the size-of-terms estimate (about
twice the largest error measured), plus 8 eps (1 + |log Psi|) for the
error of the normalising constant itself (at most 4.4 eps (1 + |log Psi|)
on scans of 1.55e5 points). Against an arbitrary-precision reference it is
at most 1.4e-10 at 1005 such points on 390 triples, and 1.1e-10 at 945
points placed from just below the gate to ten times it on 244 random
triples with alpha < 1 and gamma / sqrt(beta) from -3e4 to -20, where
0.1.1 is up to 4.0e-7 off, and 6.1e-11 at 600 points from just above the
gate to ten times it on 150 random triples with alpha from 5 to 200 and
gamma / sqrt(beta) from -40 to -1. 0.1.0 was right here.
In the same way pmhn(q, log.p = TRUE) is now formed in logarithms
throughout. Outside the special cases 0.1.1 took the logarithm of F:
near F = 1, at the 1e-10 upper quantile of (2, 1, 1), it returned
-1.00000008e-10 for -1e-10, and wherever F fell below the smallest
positive double, about 4.9e-324, it returned -Inf:
pmhn(75, 2, 1, 300, log.p = TRUE) was -Inf and is -5631.28. Where F
is a subnormal double that logarithm kept only part of its digits.
For a shape near 1e-300 at a small beta, that lower tail needed one more
change. Next to the origin the integral of the density is about
1 / alpha, and the lower-tail quadrature measured it against the kernel
at q, which carries the factor 1 / q: at
pmhn(1e145, 1e-300, 1e-300, 1001e-150, log.p = TRUE) their ratio is
exp(1024), past the largest double, so that the integral would
overflow. Where that part exceeds the kernel by more than exp(600) it is
now the measure instead, which leaves the integral as it is; the call
returns -249803.83 (0.1.1: -Inf).
For a small shape pmhn(q, log.p = TRUE) lost the relative accuracy of
the lower tail below the peak, where F is close to 1: its logarithm,
about -P(X > q), came from terms of the size of log Psi,
(alpha / 2) log(beta) and (alpha - 1) log(q) that cancel. Under a
small tilt the error passed 1e-9 below a shape of about 1e-6 (7.5e-11 at
1e-5, 1.1e-9 at 1e-6), and under a large tilt below a smaller shape.
pmhn(2.5, 1e-100, 1, 5, log.p = TRUE) was 0 where it is
-1.52e-98, pmhn(2.5, 1e-12, 1, 5, log.p = TRUE) 1.1e-4 off, and
pmhn(2.67e-308, 2^-1022, 1, -0.5, log.p = TRUE) -1.1e-13 where it is
-1.57e-305. Wherever that error could pass 1e-9 of the logarithm and the
error of the upper tail, carried into it, is estimated at less than half
as much, it is now log(1 - P(X > q)) from the upper tail, as the
logarithm of the upper tail is formed from the lower one on the other side
(see above). On scans with shapes from 2.2e-308 to 1e-4, tilts from -1e5
to 45 and beta from 1e-300 to 1e300 it is within 1e-11 of independent
quadratures where it is taken (see the next entry for a peak below the
smallest normal double). Where it is taken,
the upper tail is evaluated as well as the lower one: 5.0 microseconds a
point against 3.0 at (1e-10, 1, 0.5), and 4.0 against 2.6 at
(1e-100, 1, 5).
Where the peak of the density of log X lies below the smallest normal
double, 2.2e-308, as it can for a tiny shape under a strong negative tilt
at a large beta, the upper-tail quadrature took the slope of its
integrand at the peak from the peak's stationarity, which holds only at
the exact peak, while a subnormal peak is resolved only to the spacing of
the doubles there: pmhn(1.482e-323, 1.41e-217, 7.65e208, -37.47 * sqrt(7.65e208), lower.tail = FALSE, log.p = TRUE) was 3.5e-7
off. The slope is now formed directly there. On 200 such points the
logarithm of the upper tail is within 2.2e-16 of an independent
quadrature, and that of the lower tail within 1e-13.
qmhn() solved a lower-tail p on F(x) = p unless 1 - p was below
1e-7, while F near 1 is resolved only to the rounding of 1 and, for a
small shape, to the error of the lower tail above: a relative error of
that divided by 1 - p. qmhn(log1p(-1e-6), 1e-8, 1, 0.5, log.p = TRUE)
was 1.3e-7 of itself off, and at (2, 1, 0.5) the quantile of
p = 1 - 1e-7 left P(X > q) 2.1e-9 of 1 - p away from it, 1 - p
taken as the double it is. Wherever that error could pass 1e-9 of
1 - p and the upper tail resolves it more accurately, the quantile is
now the root of P(X > x) = 1 - p, with 1 - p formed exactly: within
about 1e-5 of 1 at parameters of order one, and farther for a small or
large shape, a large tilt or an extreme beta.
At (2, 1, 0.5) that difference is now 4e-15.
qmhn(1e-100, 2, 1, 3) returned 1.62e-40, whose probability is 5.2e-82; it
now returns 7.11e-50.
qmhn(p, lower.tail = FALSE) for p close to 1 solved on the logarithm
of the upper tail, which is about -(1 - p) there and carries the
absolute error described above. That is a relative error in 1 - p, and
at 1 - p = 1e-6 it moved the quantile by 8.2e-5 of a standard deviation
at (9e5, 1, -500), 7.3e-5 at (2e4, 1e-300, 5e-151), 4.6e-5 at
(0.3, 1, 300) and 7.4e-4 at (0.5, 1, 900). At the last two most of
that error is the error of the normalising constant under a strong
positive tilt (2.3e-10 and 3.6e-9 in its logarithm, fixed above), which
this equation turns into a relative error in 1 - p; with the constant
fixed, the equation alone leaves 8.9e-7 and 3.4e-6. Where that relative
error could pass 1e-9, the quantile is now the root of F(x) = 1 - p,
found directly where F comes from a quadrature, and by Newton steps from
the upper-tail root where it comes from the series of Sun et al. (within
1e-7 of 1, where needed, from the root for a 1 - p the upper tail still
resolves). On 141 triples at upper-tail p from 0.7 to 1 - 1e-12 it is
within 1.3e-9 of a standard deviation of an arbitrary-precision reference,
where 0.1.1 is up to 8.8e-3 off down to 1 - 1e-7 (at
(0.5, 1e-300, 1e-147), again mostly the error of the normalising
constant; the equation alone left up to 2.6e-3, at (2e4, 1e300, 5e151)).
Its relative error is that of 1 - p divided by x f(x) / F(x), about
alpha next to the origin: on 178 other triples it is up to 6.8e-10, at
alpha of about 0.015, against up to 2.9e-5 in 0.1.1 down to 1 - 1e-7;
closer to 1, 0.1.1 missed by 2.7e-3 or more at a tenth of the points.
Where the quantile lies below the smallest positive double, as it can for
alpha < 1 under a negative tilt, 0.1.1 returned a positive value
(qmhn(log1p(-1e-35), 0.0677, 1, -267.08, lower.tail = FALSE, log.p = TRUE) was 3.9e-229); it now returns 0, as the lower-tail
quantile of 1 - p does.
At gamma = 0 the square-root-Gamma forms squared x, which underflows
below 1.5e-154: dmhn(1e-170, 1.5, 1, 0) was Inf,
pmhn(1e-200, 0.01, 1, 0) was 0 where it is 0.0100, and
qmhn(0.01, 0.01, 1, 0) was 0 where it is 7.5e-201.
At alpha = 1 the lower tail and the quantile cancelled:
pmhn(1e-16, 1, 1, -3) was 0 where it is 3.5e-16, and
qmhn(1e-17, 1, 1, -3) was 0 where it is 2.85e-18.
When qmhn() cannot bracket a root it now returns NA, as its warning
says; 0.1.1 returned NaN. A qmhn() call with scalar parameters can now
be interrupted.
For a shape below about 1e-11 under a positive tilt, the upper tail at a
very small q left out the stretch next to the origin, which holds most
of the tail there: pmhn(1e-200, 1e-20, 1, 1, lower.tail = FALSE) returned
3.29e-19 where the value is 4.61e-18, and
pmhn(1e-30, 1e-16, 1, 1, lower.tail = FALSE) 6.954e-15 where it is
7.003e-15. qmhn() with lower.tail = FALSE inherited the error:
qmhn(1e-12, 1e-14, 1, 1, lower.tail = FALSE) returned 9.660e-44 and now
returns 9.616e-44.
RTDR drew 0 too often at a shape just above the smallest normal double
where the law puts a fair share of its mass below the smallest positive
double: at alpha = 2^-1022 and a tilt of 53.333, where that share is
1/2, rmhn() drew 0 for 50.8 percent of the draws (15 standard errors
over 1e6 draws), and at 1.25 * 2^-1022 and 53.325, where it is 1/2
too, for 50.3 percent. On the log
axis the leftmost piece of the envelope falls like exp(alpha y), and a
draw from it, a + log(u) / alpha, passed the largest double and became
-Inf once |log u| exceeded alpha times the largest double, which at
2^-1022 happens for u below e^-4; the acceptance test then compared
-Inf with -Inf and accepted every such draw. Such a draw is now tested
at the largest negative double, where the ratio of the law to the
envelope has already reached the value it keeps below y = -745, and is
accepted with that probability, so that the acceptance rate is again the
ratio of the two areas. Under R's default generator, whose uniform
variates are not below 2^-33, this can happen only below a shape of about
1.3e-307; under another generator it can reach somewhat larger shapes,
and is handled the same way.
Sun et al. Algorithm 3 searched for its matching point with an absolute
stopping width, so once beta passed about 1e24 its acceptance rate
collapsed and draws came back as NaN: rmhn(2000, 50, 1e30, -3e15, method = "sun") returned NaN for all but a handful of its draws, with a
warning. The default method used Algorithm 3 for many of these laws, so it
was affected too: rmhn(500, 10, 1e30, -5e15) returned NaN for about 70%
of its draws. The search now runs at beta = 1 and is rescaled.
Algorithm 1 drew from too wide a law at a very large shape: its standard
deviation was about 22% too large at (1e20, 1, 1e10) and 46% at
(1e28, 1, 1e10), and a Kolmogorov-Smirnov test of 20,000 draws against
the normal limit gave p below 1e-40. The standard deviation of 20,000
draws is now within 1% of the limit's. For beta above about 1e154 it
returned only NaN, as in rmhn(2000, 2, 1e160, 1e80).
For small shapes a Gamma variate that underflowed became an exact zero: at
(0.01, 1, 0) 2.4% of the draws were 0, where the law puts 0.058% of its
mass below the point at which a double rounds to 0, and at alpha = 0.005
15.7% against 2.41%. Such a variate is now replaced by a draw from the exact
tail of the law.
For alpha < 1 with a strong positive tilt, RTDR built its envelope on an
absolute log axis. From gamma / sqrt(beta) of about 1e13 its draws fell
on a lattice: at (0.7, 1, 1e14) 20,000 draws took only about 30 distinct
values, and a Kolmogorov-Smirnov test against the normal limit gave p
below 1e-100. They now take some 570 values, the doubles 2^-7 apart near
the mode, with the frequencies of the normal limit rounded to them. For
1/2 <= alpha < 1 it stopped from about 3e14 with an error that called the
parameter values a defect, and for alpha < 1/2 from about 3.3e14.
All regions now build and sample their pieces in the offset from the mode,
and stop only where the law is narrower than the spacing of the doubles at
its mode (see below).
On the default path for gamma < 0, RTDR drew from the wrong law, without a
warning, once the mode of sqrt(beta) X fell below about 1e-27:
rmhn(20000, 2, 1e-56, -1) had mean 1.17 where the law, Gamma(2, 1) to
double precision, has mean 2 (Kolmogorov-Smirnov p = 0; the draws now
have mean 2.0), and rmhn(20000, 4, 1, -1e30) stopped with an error that
called itself a defect.
For alpha > 1 under a tilt gamma / sqrt(beta) near -alpha / e, where
the mode of sqrt(beta) X lies near e, RTDR drew from the wrong law
without a warning, or stopped with an error that called itself a defect,
from a shape of about 3e15: rmhn(20000, 1e17, 1, -1e17 / exp(1), method = "rtdr") had a mean about 0.75 standard deviations below the law's
(Kolmogorov-Smirnov p = 0 against Algorithm 3), and
rmhn(50, 3e17, 1, -3e17 / exp(1), method = "rtdr") stopped. Its setup
chose whether to build the envelope relative to the mode from the size of
the log-density at the mode, which there is the small difference of terms
of order alpha and carries their rounding error; it now chooses from the
size of those terms. Over 81 parameter sets with shapes from 1e14 to 5e17
and tilts within 1% of -alpha / e, 0.1.1 drew the wrong law in at least 21
and stopped in 8; RTDR's draws now agree with those of Algorithm 3 in all 81
(two-sample Kolmogorov-Smirnov tests of 20,000 draws each, smallest p
0.0013). The default method of 0.1.1 used Algorithm 3 there; the default
method now uses RTDR there at any number of draws (see Performance).
rmhn(20, 20, 1e-300, -1e300) returned only NaN.
Under method = "rtdr", a law of tiny shape under a strong negative tilt,
whose mass lies almost wholly below the smallest positive double, stopped
with an error that called itself a defect, as in
rmhn(10, 1e-17, 1, -1e300, method = "rtdr"). It is now drawn by Algorithm
3, which returns zeros there.
The help pages of dmhn(), pmhn() and qmhn() now state the range over
which their results are accurate and what is returned beyond it, and that
of rmhn() which laws double precision lets it sample and what it returns
or stops with where it cannot. For pmhn() and qmhn() that range is
alpha and |gamma| / sqrt(beta) up to 1e15 with beta in
[1e-300, 1e300]. Beyond it the law can be narrower than the spacing of
the doubles near its peak, and the probabilities are those of the exact law
at the doubles given: pmhn(mhn_mode(2, 1, 1e16), 2, 1, 1e16) is now 0.5
(0.1.1: 0.585), and pmhn(mhn_mode(2, 1e26, 1e17), 2, 1e26, 1e17) 0.5
(0.1.1: 0). Under a nonzero tilt, from a shape of about 1e22, 0.1.1's
pmhn() also returned NA or stopped, with an error from the quadrature
library, a length error or std::bad_alloc; at gamma = 0 it did not
stop.
Past |gamma| / sqrt(beta) of about 1e154, where gamma^2 overflows, the
moment helpers stopped with an error from the quadrature library, and
mhn_mode(2, 1, 1e200) was Inf. They now return the values of the limit
laws: mhn_mean(2, 1, 1e200) and mhn_var(2, 1, 1e200) are 5e199 and 0.5
(the normal limit), and mhn_mean(2, 1, -1e200) and
mhn_skewness(2, 1, -1e200) are 2e-200 and 1.414 (the Gamma limit).
At alpha = 1 the density was NaN where 2 beta overflows or the tilt is
close to the largest double: dmhn(1e-309, 1, 1, -1e308) and
dmhn(0.5, 1, 1e308, 1e308) are now 9.05e307 and 5.64e153.
RTDR stops with an error, saying why, only where the standard deviation of
the law is no larger than the spacing of the doubles at its mode: from
gamma / sqrt(beta) of about 9e15 (2^53 at beta = 1), and under
method = "rtdr" with gamma <= 0 from a shape between 2^103 and
2^106. So rmhn(3, 0.3, 1, 8e15), which 0.1.1 refused, now samples,
while rmhn(3, 1.5e31, 1, 0, method = "rtdr"), whose draws 0.1.1 centred
on the mode rounded to a double (0.035 standard deviations from the true
centre), now stops; the default method still samples it.
method = "sun" now stops with an error for a negative tilt (gamma < 0,
other than a negligible one, which is drawn as at gamma = 0) above a shape
of 2^59 (about 5.8e17), where R's rgamma() can shift the
variates it supplies by more than 1e-7 of a standard deviation; the default
method uses RTDR there. Above alpha = 2^60 a negligible tilt is drawn from
the Normal proposal of Algorithm 1 centred at the mode (by RTDR under
method = "rtdr") instead of as the square root of a Gamma variate.
With beta far outside [1e-300, 1e300], two intermediate quantities of
the closed forms for the mode and for the points the distribution function
starts from left the double range although the results do not. Where
alpha and beta both exceed about 2.2e307, sqrt(8 beta (alpha - 1))
overflows: mhn_mode(1e308, 1e308, -1e308) is now 0.5 (0.1.1: NaN), and
pmhn(0.9 * (1 + sqrt(5)) / 2, 1e308, 5e307, 1e308), a point below a law
that sits on its mode 1.618, is 0 (0.1.1 stopped with an error from the
quadrature library). Where gamma^2 and 8 beta |alpha - 1| both fall
below the smallest normal double (for alpha = 2, a beta below about
2.8e-309 with |gamma| below 1.5e-154), they lose digits:
mhn_mode(2, 1e-322, 3e-161) is now 1.799249e161 (0.1.1: 1.799002e161),
and pmhn(0.9 / sqrt(1e-320), 150, 1e-320, 1e-170, log.p = TRUE) is
-268.49 (0.1.1: NaN). On 1287 parameter triples with beta in
[1e-300, 1e300] every function and every seeded draw is unchanged by
this, except mhn_mode() at one triple (alpha = 1 + 1e-9,
beta = 1e-300), in its last digit.
dmhn(), pmhn(), qmhn() and rmhn() treat invalid arguments as the
base R distribution functions do (?dgamma: "Invalid arguments will result
in return value NaN, with a warning"), element by element:
alpha or
beta not positive and finite, or gamma not finite) is NaN, with the
warning "NaNs produced" (rmhn(): "NAs produced"). 0.1.1 stopped the
whole call, in rmhn() for a zero or negative alpha or beta.NA or NaN parameter gives NA or NaN without a warning from
dmhn(), pmhn() and qmhn(), where 0.1.1 stopped the whole call, and
NaN with the warning "NAs produced" from rmhn(), where 0.1.1 gave NA
without one. An NA input gives NA and a NaN input NaN;
dmhn() had relied on the exponential to keep an NA as NA, which R
does not promise.alpha, beta or gamma gives numeric(0), and in
rmhn() NA draws with a warning, where 0.1.1 stopped.The other elements are unaffected, and an invalid element of rmhn() draws
no random numbers. The moment functions and mhn_mode(), which take a
single triple, still stop on an invalid one.
A factor, a character string, a complex number, a list or NULL given to
dmhn(), pmhn(), qmhn() or rmhn() is an error, as in dnorm();
0.1.1 used a factor through its level codes, so that dmhn(factor("3"))
was the density at 1, and a string through as.numeric(). A Date or
another vector stored as numbers is taken as its numbers, as before. The
moment functions and mhn_mode() say which argument is not a single
number, where 0.1.1 gave Rcpp's "Expecting a single value" or "Not
compatible with requested type".
qmhn() now warns "NaNs produced" for a probability outside [0, 1], as
the base R quantile functions do.
dmhn(), pmhn() and qmhn() keep the attributes of their first argument,
such as names and dim, as the base R d/p/q functions do.
rmhn() uses only the elements of alpha, beta and gamma that its
draws use, the first n, and ignores the rest, as rnorm() does:
rmhn(1, alpha = c(2, -1)) stopped with "alpha must be positive" and now
draws from alpha = 2.
rmhn() reads n as rnorm() does: a zero-length n asks for no
variates, where 0.1.1 stopped, and NULL and a negative fraction such as
-0.5 are errors, where 0.1.1 returned numeric(0) for the latter. An
n that is NA, NaN, infinite, beyond the integer range or a string
that is not a number is reported as such rather than as negative, and one
beyond the integer range without a coercion warning on the way.
Messages from internal routines (std::bad_alloc, quadrature-library
errors) no longer reach the user, and the errors of rmhn() say what went
wrong at which parameter values instead of calling a representability limit
a defect. The help page of rmhn() documents the NaN draws and the errors
it can return.
The warnings of the compiled code (dmhn(), pmhn(), qmhn(), rmhn())
are raised through R's warning() under Rcpp's unwind protection (its
default from Rcpp 1.0.10), and carry no call. A warning that is caught by
tryCatch(..., warning = ), or turned
into an error by options(warn = 2), does not return to the compiled code;
raised by Rf_warning(), as Rcpp::warning() raises it, it would leave by
a long jump that skips the C++ destructors, and the vectors of the call
would stay allocated for the rest of the session.
mhn now requires Rcpp 1.0.13: Imports and LinkingTo ask for
Rcpp (>= 1.0.13), where 0.1.1 asked for 1.0.7 in Imports and for any
version in LinkingTo. Rcpp turns on the unwind protection above by default
from 1.0.10, and from 1.0.13 its List::create() takes more than 20
elements, as an internal routine of mhn needs; mhn does not compile against
Rcpp 1.0.12.
Measured on an Apple M4 Max against 0.1.1, as medians of repeated calls:
dmhn() at the default parameters is about 4 times as fast (1e5 values in
0.98 ms against 4.3 ms), and 1.7 times as fast with recycled parameter
vectors.
pmhn() for q above the median is 3 to 12 times as fast (5 against 17
microseconds per value at (2.5, 1, -8)), and about 160 times at
(2.5, 1, 30); for |gamma| / sqrt(beta) above 20 it is faster on both
sides of the median (0.36 against 520 microseconds per value below it at
(2, 1, -50)); and at a large shape pmhn(71, 1e4, 1, 1) takes a few
microseconds against about 8 milliseconds.
qmhn() is 1.8 times as fast at (2.5, 1, 1.5), and at alpha = 1 5.6
times with recycled parameters and about 150 times for
qmhn(p, 1, 1, -5); qmhn(0.5, 1.5, 1, 300) takes about 30 microseconds
against 0.2 s, and qmhn(p, 2, 1, 2000) about 20 microseconds per
probability against 14 ms.
mhn_mean() and mhn_var() at (2, 1, 900) take 5 microseconds together,
against 17 ms.
rmhn() with short recycled parameter vectors builds one sampler setup per
distinct triple instead of one per draw: rmhn(1e4, c(2.5, 3), 1, -2) takes
0.54 ms against 3.4 ms. A single draw with method = "sun" or
method = "rtdr" given exactly takes half as long, and a batch of draws
from Algorithm 3 takes about 13% less time.
rmhn() with the default method chooses between RTDR and Sun et al.'s
Algorithms 1 and 3 by thresholds derived from timings of the current
samplers (Apple M1 Ultra), and checked by timing the code of this release.
For gamma < 0 it uses RTDR from 25 draws per
setup for alpha < 0.1 (0.1.1: from 100), from 5 for 0.1 <= alpha < 1
(0.1.1: from 25), at any number of draws for 1 < alpha < 1000, and from
10, 25 or 50 for alpha from 1000 to 1e5, from 1e5 to 1e8 and above
(0.1.1: from 25, and never for alpha >= 10). RTDR's setup walks its
two contact points in from half and about 1.5 times the mode, which from
a shape of about 1000 costs more than Algorithm 3's setup; where it starts
them next to the mode instead, under a very strong negative tilt
(gamma / sqrt(beta) < -1e10 * sqrt(alpha)) or, at moderate tilts, at a
shape beyond about 2.4e12 to 3.1e12, near a tilt of -alpha / e beyond
about 2.3e13, and at any tilt beyond about 4.5e13, the default method uses
RTDR at any number of draws. For
gamma > 0 and alpha > 1 it uses RTDR when
gamma / sqrt(beta) <= 2 * sqrt(alpha) and at least 5, 25, 50 or 100
draws (for a tilt up to sqrt(alpha)) or 25, 50, 100 or 1000 (above it)
share a setup, for alpha below 10, from 10 to 1000, from 1000 to 1e8
and above, or 5 or 25 where RTDR's setup starts next to the mode, and
Algorithm 1 otherwise (0.1.1: Algorithm 1 always). Over a paired timing
race of the two samplers in 560 cells the chosen sampler is 0.025% slower
than the faster one on average and 6.7% at worst, against 6.4% and 69%
under the rule of 0.1.1; over a second race of 1458 cells, with shapes up
to 1e14, tilts up to 3 * sqrt(alpha) and single draws within one
vectorised call (one setup per draw), it is 0.1% slower on average and
12% at worst, and never slower in the vectorised calls, against 7.6% and
47%; in three runs of inst/benchmarks/auto_dispatch.R it is at most
0.7% slower, against 65%. The draws per setup are counted as n for
scalar parameters and, where short recycled parameter vectors keep a setup
for each position within the period after which their triples repeat (a
period of at most 256 draws, recycled at least four times), as n divided
by that period, rounded down; otherwise as 1. 0.1.1 divided n by the length of the
longest parameter vector, which with recycled parameters could count
several draws per setup where every draw rebuilt one.
Slower. The figures below are per value against 0.1.1, from three rounds in which the two versions ran side by side in persistent R sessions, alternating measure by measure, on a machine carrying other load; the ratios agree between rounds to within 10% at most triples and 20% at all.
pmhn(q, lower.tail = FALSE, log.p = TRUE) also evaluates the lower
tail where the logarithm of the upper tail would lose its relative
accuracy (see above): below F of about 1e-5 at parameters of order
one, and of about 7e-4 at (50, 1, 20). There it costs about what the
lower tail costs, in either version. For |gamma| / sqrt(beta) up to
20 and alpha up to 100 that lower tail is the series of Sun et al.,
whose length grows like gamma^2 / beta with the tilt and also with
the shape, and so does the factor. With q spread evenly in log F
from 1e-6 to 1e-4 it takes 2.2 and 2.6 times as long as 0.1.1 at
(0.5, 1, -5) and (2.5, 1, 1.5), 6.5 and 8.7 times at
(10, 1, 5) and (3, 1, 10), 11 and 15 times at (0.9, 1, 12) and
(5, 1, 13), 26 to 44 times for gamma / sqrt(beta) from 15 to 19.5
at shapes from 0.01 to 10 (47 against 1.1 microseconds at
(0.5, 1, 19.5)), and 45 to 77 times at shapes from 50 to 85 (90
against 1.2 at (85, 1, -19)); with q spread evenly in log F from
1e-6 to 0.5, 1.5 to 23 times and 18 to 33 times. At the same
|gamma| / sqrt(beta) a negative tilt costs less: 7 and 12 times at
(0.5, 1, -19.5) and (3, 1, -18) from 1e-6 to 1e-4. Where the lower
tail comes from a quadrature (alpha > 100, or
|gamma| / sqrt(beta) > 20 with alpha >= 1), both tails are
quadratures, and from 1e-6 to 1e-4 it takes 1.1 to 1.9 times as long
(1.8 against 1.0 microseconds at (400, 1, -3), 1.5 times at
(1.5, 1, 25) and 1.1 at (5, 1, -40)). For alpha < 1 with a larger
tilt, where the series would cost milliseconds and the lower tail is
taken from its quadrature, checked against the upper tail, it takes 1.0
to 1.1 times as long at (0.3, 1, -30), (0.3, 1, -300) and three
triples with alpha from 0.16 to 0.39 and gamma / sqrt(beta) from
-120 to -1993, and less time under a strong positive tilt, where the
normalising constant is now cheaper (0.9 to 3.4 against 13 to 25
microseconds at (0.3, 1, 300)). With q spread evenly in F below
the median, where the upper tail alone serves, it takes as long as
0.1.1 to within 13% at all these triples.
qmhn(p, lower.tail = FALSE) for p close to 1 now works from the
lower tail (see above). Where that is the series, which the Newton
steps that refine the root sum a few times per value, it costs more:
for 20 values of 1 - p spread evenly in their logarithm, from 1e-12
to 9e-8, from 1e-7 to 1e-5 or from 1e-5 to 7e-4, it takes 1.3 to 2.6
times as long as 0.1.1 for gamma / sqrt(beta) from 15 to 19.5 (118
against 46 microseconds at (0.5, 1, 19.5) from 1e-12 to 9e-8), 2.0 to
4.0 times at shapes from 50 to 85 (236 against 59 at (85, 1, -19)),
and 1.15 times at (3, 1, -18) from 1e-7 to 1e-5. At the other triples
measured it takes at most 7% longer, and mostly less time: 74 times
less under a strong positive tilt (92 against 6800 microseconds at
(0.5, 1, 900)) and 58 times at (0.16, 1, -1993) (16 against 964).
A single call of dmhn() with scalar arguments takes up to about 0.1
microseconds longer, 2 to 11% in the measurements, for keeping the
attributes of x (see Interface).
A shape within sqrt(.Machine$double.eps) of 1 or a tilt below it,
other than exactly 1 or 0, now costs what the general case costs next
to it: pmhn() takes 4.2 microseconds per value at
(1 + 1e-12, 1, -1) and qmhn() 40, against 0.03 for both at
alpha = 1.
For |gamma| / sqrt(beta) up to 20, the lower tail of pmhn() below
the median, which is still summed from the series, takes up to 8%
longer at some of these triples (8% at (15.9, 1, -1.69), 7% at
(50, 1, -15) and (85, 1, -19)), and up to a third less at others.
rmhn() for gamma > 0 and alpha > 1 with a parameter triple that
changes at every draw (one Algorithm 1 setup per draw, as in a vectorised
Gibbs update) takes 9 to 25% longer per draw than 0.1.1 at shapes from 2.5
to 1e6, 7 to 14 nanoseconds more, 3 to 18% longer at shapes of 1.1 and
1.5, and up to 13% longer at shapes from 1e8 to 1e12, where near a tilt
gamma / sqrt(beta) of sqrt(alpha) it takes up to 20% less time instead
(vectors of 2000 triples, 8 to 12 rounds). The setup of Algorithm 1 now
forms the constants of the acceptance test of its proposal and checks that
they stay within the double range.
Where 0.1.1 was wrong, the right answer can also cost more:
mhn_kurtosis(1000, 1, 0) takes 11 microseconds against 6.
In the two timing races of the default method's choice (see above),
which time both samplers of this version, the sampler rmhn() now
picks is at most 2.3% slower than the one the rule of 0.1.1 picks,
except at alpha = 1.5 with gamma / sqrt(beta) = 2 * sqrt(alpha),
where from 25 draws per setup it picks RTDR, 4.6 to 5.4% slower than
Algorithm 1.
Nothing else was measured to be slower beyond run-to-run noise: 20 further workloads covering the density, both tails, the quantile, the moments and the samplers (single draws, batches of 1e5 from each sampler and a call with recycled parameters) took at most 7% longer.
The help pages state what each function returns and over which range its
results are accurate, and leave the numerical methods to the theory
vignette. The main text of the vignette keeps to the theory and to how
method = "auto" chooses a sampler; a new Appendix A gives the numerical
methods, the handling of extreme parameters and the timings behind
method = "auto", and Appendix B the misprints in the source papers.
The NEWS entry of 0.1.1 attributed to Sun et al. (2023) a misprint in the
constant K_1 of the first part of Theorem 1 (printed "a."), and the
theory vignette of 0.1.1 one in the Gamma argument of the mixture weights
of Lemma 7. The printed article contains neither: that part of Theorem 1
prints the numerator of the base of K_1 as sqrt(beta) (alpha - 1) and
Lemma 7 prints Gamma(alpha/2 + i/2). The NEWS entry is corrected below,
and the vignette now lists the six misprints the article does contain in
the results it draws on that affect what is computed from them (it has
other slips there, which do not); the package uses the correct form of
each result it implements (the weights of Lemma 7 serve only
Algorithm 2, which it does not implement):
I(K_1(delta, ...) <= K_2(mu, ...)), with the first arguments of K_1
and K_2 exchanged;1 / (2 beta^(alpha/2)),
so they do not sum to one;alpha x and (alpha + 1) x inside the
discriminants of C_1 and C_2, where the derivation gives
alpha x^2 and (alpha + 1) x^2;(beta m + |gamma| alpha) / (2 beta m + |gamma|), where the display that
introduces the proposal has
alpha (beta m + |gamma|) / (2 beta m + |gamma|);l'' of Supplementary Section 2.17 prints its first term
without the minus sign that differentiating l' gives.The entry for 0.1.1 below is corrected where it did not describe 0.1.1 or
0.1.0 as released: the counts of values and draws that changed, the scope
of the RTDR overflow, the acceptance rate and the behaviour of RTDR at a
steep tilt, four failures it attributed to 0.1.0 that 0.1.0 did not show,
the cost of the series at gamma = -1e4, the value qmhn() returns when
it cannot bracket a root, the loops that check for an interrupt, and the
size of the test suite.
The help page of rmhn() said that a -Inf in alpha or beta gives an
NA draw, whereas 0.1.1 stopped the call for a zero or negative alpha
or beta, -Inf included; 0.1.2 gives NaN with a warning (see
Interface). The help page now also says how the number of variates per
setup that steers the default method is counted, and that the switch to
RTDR above alpha = 2^59 concerns a negative tilt only.
The help pages describe beta as the rate-type coefficient of
exp(-beta x^2) and gamma as the coefficient of the tilt exp(gamma x)
(0.1.1 called them a scale and a location parameter), and the theory
vignette describes the normalising constant, the distribution function and
the quantile as they are now computed.
Appendix A of the theory vignette gives the accuracy of the ordinary forms
of pmhn() as measured: below 1e-8 within 10 standard deviations of the
peak, and about 3e-8 at 36 standard deviations in the upper tail close to
the bound at which the forms measured from the peak take over. The help
page of pmhn() names a third case
outside the stated range, beta from about 2e307, in
which a point far above the peak gives NA with a warning or -Inf.
The help page of qmhn() says that a lower-tail p is solved on the
series where pmhn() uses it, which for alpha < 1 with |z| from about
100 can take from a tenth of a second to seconds for one quantile, and that
1 - p with lower.tail = FALSE gives the same quantile at once. It says
what is returned where the log upper tail of pmhn() is -Inf for beta
from about 2e307.
The help pages of dmhn(), pmhn() and qmhn() say when the normalising
constant is reused across a vector of parameters: for each position of a
recycling cycle of at most 256, and otherwise for each run of consecutive
elements with the same parameters. The help page of dmhn() also says
that log = TRUE returns -Inf in place of a finite log density for a
positive tilt whose mode lies beyond the largest double.
The help pages of mhn_skewness() and mhn_kurtosis() state the accuracy
measured against independent quadratures (an absolute error below 1e-9 for
shapes from 1 to 60 and tilts from -20 to 12, the largest, 3e-10 and 6e-10,
for a shape between 2 and 3 under a positive tilt), and how a subnormal
shape is served.
citation("mhn") gives the package's DOI, 10.32614/CRAN.package.mhn, and
the reference to Gao & Wang (2025) in the help pages, the README and the
vignettes says it is an advance online publication, pages 1-25.
R CMD check, where 26 heavier blocks, two timing
tests, a memory test and a comparison to the last bit are skipped.The result files that 0.1.1 installed under benchmarks/results and
audits/results are no longer part of the package:
system.file("benchmarks", "results", package = "mhn") and
system.file("audits", "results", package = "mhn") are now "". The
scripts under benchmarks/ and audits/ and the example under
examples/ remain.
The benchmark and audit scripts and the example write their results only
to the directory named by their output-directory variable (MHN_BENCH_OUTDIR,
MHN_GW_OUTDIR, MHN_GOF_OUTDIR, MHN_VMF_OUTDIR), and otherwise only
print a summary, where 0.1.1 created a results directory inside the
working directory. Next to the results they record the mhn version and
source commit, R, the platform, CPU, compiler flags, BLAS and the run times
(benchmarks/provenance.R). auto_dispatch.R takes its seed from
MHN_BENCH_SEED; replicate_gw_tables.R prints the ratio of the two
samplers' times beside Gao & Wang's ratio, their times not being
comparable with the package's; and vmf_gibbs.R no longer sums every
autocorrelation into the effective sample size when the one at lag 1 is
negative.
series_breakdown_audit.R forms the incomplete gamma values of its
500-bit reference from their logarithms. 0.1.1 took them in double
precision, where under a strong positive tilt they underflow to zero, so
that the audit reported breakdowns of the double series that were errors
of the reference. It also measures the error of a logarithm as the
absolute difference, which is the relative error of the value.
Bug-fix and maintenance release.
This release corrects a number of numerical defects reachable from the public API. Most returned a wrong number silently rather than raising an error, and most require parameter values outside the ranges 0.1.0 was tested over. Each fix below carries a regression test at the values that failed.
Numbers change wherever a defect is fixed, and that includes part of the
ordinary parameter range --- not only the extremes. Everything below was
measured against 0.1.0 as published, over a grid of 330 parameter triples
spanning alpha in [0.05, 20], beta in [0.01, 100] and gamma in
[-1e6, 1e6].
Deterministic functions. Over the 130 grid cells with
|gamma| / sqrt(beta) at most 10, dmhn(), pmhn(), mhn_mean() and
mhn_var() agree with 0.1.0 to twelve significant digits in 100 cells and
differ in 30. The largest change is 8.2e-6 relative, in mhn_var(); dmhn()
and pmhn() move by at most 3.1e-9 and mhn_mean() by 1.6e-7. In every case
checked against an independent arbitrary-precision reference the 0.1.1 value
is the accurate one --- the Lemma 2c variance already cancels measurably at
|gamma| / sqrt(beta) = 10.
Outside that range the differences are far larger, and are the defects listed
below.
rmhn() is not seed-compatible with 0.1.0. For a fixed seed, 163 of the
330 cells return a different sequence of draws; 51 of those are in the ordinary
range above. Draws at gamma = 0 are unchanged, since that case is a
closed-form transform of rgamma(). There are two reasons a sequence changes,
and they matter differently:
(alpha, beta, gamma) = (0.3, 1, 100) a
Kolmogorov-Smirnov test of the 0.1.0 sample against the true distribution
gives p = 0; the same test on 0.1.1 gives p = 0.51.rmhn() is a rejection sampler: correcting a constant in the envelope changes
which variates are proposed and which are accepted, so the stream is consumed
differently. At (10, 1, 1) every one of 200,000 draws from a fixed seed
moves, and yet a Kolmogorov-Smirnov test against pmhn() gives p = 0.78
on the 0.1.0 sample and p = 0.70 on the 0.1.1 one, with both sample means
within 0.05% of the true 2.4387. Analyses in this category remain valid;
only bit-for-bit reproduction of a stored sample is lost.rmhn(method = "rtdr"), and hence the default method = "auto" where it
routes there, drew biased samples for alpha < 1 and gamma > 0 outside
the corner alpha < 1/2, gamma / sqrt(beta) > 2 * (1 - sqrt(1 - 2 * alpha)):
on the log axis the density is only T_{-1/2}-concave in that region, not
log-concave, so the previous log-tangent envelope did not dominate it and
over-weighted large values. The envelope now follows Gao & Wang (2025,
Section 3.2 and Appendix B), using a T_{-1/2} (inverse-square) tangent hat
when gamma > 0 and the log-tangent hat only when gamma <= 0; the
inflection-point envelope used in that corner was not affected.
RTDR drew badly biased samples wherever an exponential envelope piece rose
by more than t = 709 in log scale, because the piece was evaluated by
forming exp(t) before dividing it out. For alpha < 1/2 this is region
D's left-tangent piece once gamma / sqrt(beta) exceeds about 76
(t = 1068 at 100), under both the default method and method = "rtdr".
For alpha > 1 it is region A's left-tail piece under method = "rtdr"
once gamma / sqrt(beta) exceeds about 530. Once alpha exceeds about 5e4
it also happens at small and negative tilts, which the default method
reached for gamma < 0 whenever 25 or more draws shared one parameter set.
The piece area became infinite, so piece selection always fell through to
the last piece, and the widest piece proposed nothing inside its own
support. At (0.3, 1, 100) a Kolmogorov-Smirnov test against pmhn() gave
p = 1.2e-121. Both the area and the inverse CDF now factor the dominant
exponential out instead of forming it.
The RTDR envelope was decided by rounding once the standardised tilt
gamma / sqrt(beta) passed about 1e8, and answered differently on different
machines: at (0.7, 1e-8, 1e8) macOS returned 50 NaN while Linux and
Windows raised an error, from the same input, and the outcome was not even
monotone in gamma. Near the mode the density carries the constant
gamma_norm^2 / 4 --- 2.5e23 at that tilt --- while every quantity the
envelope reads out of it is an O(1) difference of two such values, so the
contact drop was quantised to multiples of 3.4e7 against a target of
log 4. All three envelope regions carried it. The ordinate is now measured
from the mode, past a conditioning threshold, through an identity that only
ever adds same-signed terms, and the accept/reject test is measured the same
way --- without that second half the sampler draws from its own hat.
Fixing the arithmetic exposed three further defects that it does not itself
cure. The contact searches started a fixed distance from the mode and needed
more Newton steps than their budget allowed --- 31 against 30 at
gamma / sqrt(beta) = 1e10 --- so they returned a contact point whose drop
was hundreds instead of one; they now start where Gao & Wang (2025, Sec. 2)
recommend starting their Eq. (8) iteration. The secant count ceil(rho) was
narrowed into an int, which is undefined behaviour past a tilt of about
1.8e9. And the plateau piece was dropped without complaint whenever its
edges crossed, which is how the sampler's support came to exclude the mode.
This fix leaves draws below gamma / sqrt(beta) = 5e6 unchanged. Above it
they can change, and the new ones are the correct ones: against the
distribution's exact Gaussian limit the Kolmogorov-Smirnov statistic is
0.004 where it had reached 1.000, and the sample standard deviation is
within 0.5% of the truth where it had been 1593 times too large. The
acceptance rate --- which Gao & Wang bound below by 1/e and which no
goodness-of-fit test can see, because an envelope that is merely too wide
still samples the right law --- now holds steady from a tilt of 1e6 to 1e13:
between 0.82 and 0.89 for alpha > 1 and for alpha < 1/2, where it had
fallen to 0.001, and about 0.76 for 1/2 <= alpha < 1 with a positive tilt,
the rate of the inverse-square hat described above.
Once the law is narrower than the spacing of the doubles at its mode, no
sampler working in these coordinates can proceed, and RTDR stops with an
error instead of returning draws that are mostly NaN. Short of that point
this release does not yet sample every law correctly: for alpha < 1 the
draws fall on a coarse lattice from a tilt of about 1e13 (27 distinct values
in 2000 draws at (0.7, 1, 1e14)) until the sampler stops, near 3.3e14 for
alpha < 1/2, and for 1/2 <= alpha < 1 the error raised from about 3e14
calls the parameter values a defect.
rmhn(200, 1.5, 0.01, 3162.3) returned 200 NaN, behind 200 warnings that
the Algorithm 1 sampler had exhausted its retry budget. Sun et al. (2023)
require the Algorithm 1 proposal rate to lie in (0, beta) (Theorem 1,
first part printed "a."), and the delta_opt printed in its second part
"a." does so for every gamma > 0. Evaluated as printed in double
precision, however, it subtracts two quantities of size gamma^2, and for
a large tilt the rounded result leaves the interval --- reaching 0, then
going negative, then snapping to beta. A negative rate made the branch
test NaN < 0 select the sqrt-Gamma proposal, and rgamma was handed a
negative scale. Rationalising the expression removes every subtraction.
At alpha = 1 with gamma = -1e8 the sampler returned 81572 negative values
out of 2e5 --- outside the (0, Inf) support --- and collapsed the sample
onto 13 distinct doubles. Robert (1995) proposes y = a + e with
a = -mu/sigma, and the code returned mu + sigma * y; since
sigma * a = -mu that is sigma * e with the mean removed analytically,
but written literally it cancels. It now returns sigma * e.
The RTDR region A left-contact search started at max(m/2, 1e-6). Whenever
the mode was below 2e-6 that floor put the start on the wrong side of it,
and rmhn(200, 4, 1, -1e7) failed outright on the default path.
The RTDR right-contact searches ran on the log axis, where the density
carries exp(2y) and overflows past y = 354. The starting offset Gao &
Wang (2025, Sec. 2) recommend for their Eq. (8) iteration lands beyond that
when the mode is flat --- at alpha = 1e-5, 447 above the mode under a
negative tilt and 372 under a small positive one --- and in region D, which
started one unit above the mode, a single Newton step could land there. The
search returned a non-finite slope, which the sign guard admitted because
NaN fails every comparison, and the plateau piece was dropped with no
error: the sampler's support no longer contained the mode or the right
tail. At (3e-4, 1, 0.01) it drew nothing above 0.005, where the true
probability is 0.0015.
Corrected the Algorithm 1 envelope constant K_1. Version 0.1.0 had
implemented the numerator of its power base as sqrt(beta (alpha-1)),
where the first part of Sun et al. (2023) Theorem 1 (printed "a.") prints
sqrt(beta) (alpha-1). The sampled law was always correct --- both
proposals are valid envelopes and each acceptance test is derived
independently --- but the choice between them was not: log K_1
came out short by (alpha-1) log(alpha-1) / 2, so for alpha > 2 the
Normal proposal was selected across most of the small-tilt half of the
region. Over alpha in {4, 10} and gamma in {0.05, 0.5, 2, 20}, with
beta = 1, the worst acceptance probability was 0.69, below the 0.8 that
Theorem 2e guarantees for alpha >= 4; it is now 0.82.
Sun Algorithm 1 formed K_1 - K_2 as a difference of two large terms;
cancelling them analytically leaves an expression with no subtraction of
like-sized quantities.
RTDR region A needs an interior mode in x, and at alpha = 1 exactly there
is none, so the envelope degenerated to a plateau of height -Inf spanning
[0, Inf] and every proposal was rejected. alpha = 1 now uses the region
whose mode is strictly positive for every alpha > 0.
The matching point Sun et al. (2023, Section 4.2) print for alpha <= 1.1,
alpha^2 / (1 + alpha), contains no beta; used as printed it was not
scale-equivariant, and the acceptance rate of Sun Algorithm 3 drifted with
the scale. It is now read as the value at beta = 1 and divided by
sqrt(beta). When the inflection heuristic left the support, the setup
discarded the inflection point entirely rather than using it, which is the
limit of Sun's own expression.
rmhn(method = "sun") refused parameter values the sun path answers in closed
form, because the pre-scan tested the raw alpha and gamma while the
dispatcher routes on the scale-free predicates.
dmhn() and pmhn() returned 0 everywhere for alpha near 2 with a large
tilt --- the two-dimensional projected-normal case (alpha = p = 2) of
Sun et al. (2023, Section 2.2).
Psi[1, z] formed exp(z^2/4) in real space, which overflows at
z = gamma/sqrt(beta) > 53.28, so the normalising constant became infinite.
The term is now built as a logarithm.
The normalising constant for gamma < 0 degraded from |z| ~ 100 and
collapsed to -Inf from |z| ~ 600, taking everything downstream with it:
dmhn() returned Inf, pmhn() returned 1 everywhere, mhn_mean() NaN.
The quadrature's upper limit extended to |z|/2, the peak of the same
expression with the opposite sign, so the mass fell inside the first fraction
of a percent of the interval. The limit now follows the decay the integrand
actually has, and the range is split at the peak.
For gamma > 0 the normalising constant had only the Lemma 10 series, which
raised an error once its truncation length passed a ceiling --- on ordinary
parameter sets, since that length grows like z^2 and (2.5, 1e-4, 20) is
already past it. The quadrature that serves gamma < 0 now serves this too.
Where the closed-form mean would difference two log Psi values of size
z^2/4, it now comes from quadrature instead.
The CDF quadrature integrated [0, x] in one piece, but for a strong tilt the
mass lies within a unit or so of the peak, so a fixed-order rule returned 0 at
the distribution's own mode. In 0.1.0 that path was unreachable, because the
series raised an error first; it became reachable once the series was given
the quadrature fallback described above, and both are fixed here. The range is
now anchored at the peak, split there, and trimmed at both ends.
pmhn(..., lower.tail = FALSE) returned exactly 0, and log.p = TRUE
returned -Inf, as soon as the survival probability fell below about 1e-16,
because it was derived as 1 - F. The upper tail is now computed directly.
pmhn() at alpha = 1 always computed the upper tail and derived the lower
as 1 - that, so pmhn(3, 1, 1, 18) gave 0 where the true value is
1.076e-17. Whichever tail is smaller is now computed directly.
qmhn(1e-4, 0.1, 1, 1) returned 1.4e-17, whose actual probability is 0.019
--- 23 orders of magnitude above the true 2.3e-40, silently. For alpha < 1
the quantile of a small p sits near p^(1/alpha), and the lower bracket
bottomed out at 1.4e-17 after 30 halvings; the solver then returned that
endpoint as though it were the root. The search now runs to the denormal
floor.
qmhn() returned +Inf for every p at alpha = 1 with
gamma <~ -11.7 sqrt(beta), where the truncation mass was formed as
1 - Phi(-mu/sigma) in linear space and underflowed to exactly 0.
mhn_kurtosis() returned -701.6 at (3, 1, 1000), for a quantity that
cannot fall below -2, and mhn_skewness() had lost its sign by
gamma = 200. Both built central moments by differencing raw ones, which
cancels once the standard deviation is small next to the mean. They are now
integrated directly, in the variable centred on the peak so that no step
forms a difference larger than the answer.
log Psi[1/2, z] carried an absolute error of 0.73 at z = -1e8, in a
quantity whose true value is -17.73: z^2/4 and log Phi(z/sqrt 2) cancel.
Substituting the scaled complementary error function cancels the exponential
analytically.
The mode formula cancels for gamma < 0; Gao & Wang (2025, Eq. 9) print the
conjugate form alongside it, and their reference code selects between the two
on the sign. mhn_mode(1.1, 1, -1e8) returned 0 instead of 1e-9. Five sites
used only the unstable form.
Both places that asked whether the tilt is negligible compared |gamma|
against sqrt(.Machine$double.eps) on an absolute scale, but the shape
depends on gamma only through Delta = gamma / sqrt(beta), since
sqrt(beta) X is MHN(alpha, 1, Delta). At beta = 4e-18 a gamma of 1e-8
is an ordinary Delta = 5, and every function returned the untilted
sqrt-Gamma answer: mhn_mean(1, 4e-18, 1e-8) gave 2.82e8 against a true
1.2503e9. Both tests now use Delta.
pmhn() with gamma = -1e4 requested about 1.4e8 series terms, taking over
ten seconds and more than a gigabyte of memory. The truncation length grows
like z^2 and had no ceiling; past one the series now returns the sentinel
that already routes the caller to quadrature.
The normalising constant for gamma < 0 was also wrong by a factor that
grew without bound as alpha fell --- 1545 at alpha = 1e-6, still 0.17%
at alpha = 0.01. A fraction 1/e of the mass of x^(alpha-1) on (0, 1)
lies below x = exp(-1/alpha), which for small alpha is below any
representable double, so no quadrature can reach it. The first stretch of
the range is now integrated in closed form. The distribution function and
the central-moment quadrature had the same gap; it was invisible while the
normalising constant shared it, and appeared as soon as that was fixed.
pmhn(..., lower.tail = FALSE) returned exactly 0 for every q at small
alpha with a strong negative tilt, and qmhn inverted that by returning its
own initial bracket --- the same number for every p, with no warning. The
quadrature range was set from the curvature at the peak, which is the wrong
scale when the range is anchored anywhere else: at beta = 1e-4 with
gamma = -1000 the Gaussian width is 894 while the kernel dies within 0.08.
The range now follows the slope at the anchor. qmhn no longer returns an
unjustified value: it warns, and returns NaN, when it cannot bracket the
root.
mhn_var() multiplies any error in the mean by about gamma^2/(2 beta): the
Lemma 2c form is exact algebra, but its leading terms are each about
alpha/(2 beta) while the answer is about alpha/gamma^2. With the degraded
normalising constant it returned 3.6e-5 at (1.5, 1, -300) against a true
1.7e-5, and beyond gamma / sqrt(beta) of about -600 it returned 0 for
alpha > 1 and alpha/(2 beta) for alpha <= 1. It is now used only
where the measured cancellation is small enough, and the second central moment
comes from quadrature elsewhere --- the closed form remains the more accurate
of the two wherever it is well conditioned.
pmhn() implemented its special-case dispatch twice and the two disagreed.
The vector-parameter path formed the truncated-normal upper tail as
1 - (1 - upper), which is exactly 0 once the upper tail falls below the
rounding of 1: pmhn(3, 1, 4, -5, lower.tail = FALSE) was 1.49e-23
element-wise and 0 vectorised.
dmhn() and pmhn() mapped a NaN input to NA, where qmhn() and the
base R d/p/q family return NaN. All three now distinguish the two on
every path.
rmhn() documents an NA draw for a non-finite parameter but guarded only
gamma; an infinite alpha or beta reached the sampler.
rmhn() took n[1] where every base R r* function reads a vector n as
length(n), so rmhn(c(10, 20, 30)) drew 10 variates instead of 3.
pmhn() and qmhn() tested lower.tail and log.p with isTRUE(), which
is stricter than base R: pmhn(1, 2, 1, 1, lower.tail = 1) returned the upper
tail, 0.5661, where lower.tail = TRUE gives 0.4339 --- and base R's
pnorm(1, lower.tail = 1) agrees with lower.tail = TRUE. All three
functions (dmhn(), pmhn(), qmhn()) now read the flag with
as.logical(), which agrees with base R for TRUE, FALSE, 1 and 0,
and reject one that resolves to NA.
Errors raised when the sampler cannot build an envelope named private C++ routines and the sign of an intermediate slope. They now report what failed, at which parameter values, and that it should be reported.
Twenty of the twenty-five compiled entry points never touch the RNG and were
still generated with an Rcpp::RNGScope, paying GetRNGstate/PutRNGstate
on every call to dmhn, pmhn, qmhn and the moment helpers. Draws are
unchanged.
qmhn() evaluated the normalising constant a second time for a triple whose
CDF state had just computed it.
rmhn() skips match.arg() when method was not supplied; it accounted for
about 55% of the cost of the single-variate call a Gibbs sweep makes.
Piece selection in the RTDR sampler rebuilt a cumulative-area table on every proposal. The table is now built once per envelope and searched by bisection.
Removed development scaffolding from the installed source: references to files not in the tarball, plan markers describing a half-finished file, and stale claims of mirroring an R implementation that no longer exists. Removed three unused functions and one unused struct field.
Added src/mhn_stable.h, which holds the better-conditioned algebraic forms
in one place rather than repeating them at the call sites that need them.
Added interrupt checking to the driver loops of dmhn(), pmhn() and
rmhn(), to the loop of qmhn() over recycled parameters, and to the two
rejection loops in the truncated-normal sampler; a long qmhn() call with
recycled parameters was uninterruptible.
Moved the retry-budget warnings out of the sampler frames. Rcpp::warning is
a bare Rf_warning, which under options(warn = 2) leaves the frame by a
long jump without running destructors while the cached envelope holds five
vectors.
R CMD check, where 40 heavier blocks are skipped,
with a regression block for each defect above at the parameter values that
failed. New coverage: non-unit beta throughout the sampler grids and a
direct check of scale equivariance; the contact-point search in isolation;
the NA/NaN contract against base R; and the accuracy audit compared
against a reference that shares no code with the package.Description field, following CRAN
reviewer feedback on unexpanded acronyms.method = "auto" dispatch for gamma < 0. Previously any
batch of 25 or more variates per setup used RTDR; benchmarking across
three independent runs showed that for alpha >= 10 the Sun et al.
(2023) Algorithm 3 has the lower per-proposal cost and wins in the batch
regime too, so auto now keeps RTDR for large batches only when
alpha < 10. This lowers the worst-case slowdown of auto relative to
the per-cell optimum from about 11% to about 4% while leaving the
common Gibbs (single-variate) path unchanged. inst/benchmarks/auto_dispatch.R
gained an alpha < 1 grid, a setup/per-proposal cost decomposition, and
a comparison of the shipped rule against the measured optimum.gamma < 0 batch cutoff from 25 to 100 variates per setup when
alpha < 0.1. The crossover between Algorithm 3 and RTDR moves to larger
batches as the shape shrinks --- it sits near 25 for alpha around 0.8 but
near 100 by alpha = 0.01 --- so the old cutoff sent 25 to 99 variates per
setup to RTDR while Algorithm 3 was still up to 10% faster there.inst/benchmarks/auto_dispatch.R
that inflated the reported median_us / iqr_us times by a factor of
about 1e6. The method = "auto" dispatch decisions are ratio-based
and were unaffected, as is rmhn(method = "auto") itself.inst/benchmarks/rmhn_gof.R, that
writes Kolmogorov-Smirnov statistics and sample-vs-theory moment
summaries across the parameter grid to a CSV._diagnostics_<date>.csv with
sessionInfo(), hardware, and mhn version provenance, matching the
audit scripts.inst/examples/vmf_gibbs.R, a self-contained Gibbs sampler for the
von Mises-Fisher concentration parameter whose full conditional is an MHN
law. It is run with
source(system.file("examples", "vmf_gibbs.R", package = "mhn")), takes
its sample size, chain length and true concentration from MHN_VMF_*
environment variables, and reports whether the 95% credible interval
covers the true concentration, and the effective sample size.Initial release.
dmhn(), pmhn(), qmhn(), and rmhn() provide density,
distribution, quantile, and random generation for the Modified
Half-Normal (MHN) distribution of Sun, Kong & Pal (2023).alpha, beta, gamma, following standard R
recycling rules.alpha, beta, gamma) triple,
so grouped inputs are evaluated significantly faster than calling the
functions inside an R loop.rmhn(..., method = "auto") (default) routes each parameter triple
to the cheapest provably-correct sampler: closed-form shortcuts for
the special cases, Sun et al. (2023, Algorithms 1 and 3) where they
win,
and the Gao & Wang (2025) Relaxed Transformed Density Rejection
(RTDR) sampler elsewhere.method = "rtdr" forces RTDR with its uniform 1/e acceptance bound.method = "sun" forces Sun Algorithm 1 (gamma > 0, alpha > 1) or
Algorithm 3 (gamma <= 0); Sun Algorithm 2 is intentionally not
implemented and an unsupported combination triggers a clear error.pmhn() uses the Sun et al. (2023, Lemma 1b) series in log space,
truncated at the Sun et al. (2023, Supplementary Lemma 10(d))
constructive bound K = max(K1, K2); the truncation residual is
bounded by the user's tolerance divided by Psi.gamma < 0 the series uses sign-separated log-sum-exp + log-
diff-exp accumulation and a runtime cancellation guard derived from
the floor set by double precision: when the relative cancellation
loss would exceed the user's tolerance, pmhn() falls back to a
peak-normalised Boost.Math quadrature (Gauss-Kronrod for alpha >= 1,
tanh-sinh for alpha < 1) of the unnormalised density.qmhn() inverts pmhn() via boost::math::tools::toms748_solve on
the bracket [sqrt(eps), E(X) + 8 sqrt(Var(X))], doubling the upper
end as needed.mhn_mean(), mhn_var(), mhn_skewness(), mhn_kurtosis(), and
mhn_mode() evaluate the closed-form / recurrence-based expressions
from Sun et al. (2023, Lemmas 2 and 3).vignette("introduction", package = "mhn") walks through every
exported function with runnable examples.vignette("theory", package = "mhn") is the theoretical companion:
it covers the MHN family and its special cases, the Fox--Wright Psi
normalising constant, Algorithms 1 and 3 of Sun et al. (2023),
the four-region
Gao & Wang RTDR construction, and the rmhn(method = "auto")
decision tree.citation("mhn") returns four bibentry objects: the package, the
Sun et al. (2023) paper, the Gao & Wang (2025) paper, and Robert (1995)
for the truncated-normal sampler used at alpha = 1.