mhn

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.

What changes relative to 0.1.1

**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:

• Unchanged: RTDR below a standardised tilt 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.
• The same sequence, with each draw changed only slightly: Sun et al. Algorithm 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.
• A different sequence: RTDR, and so the default method under a positive tilt, at a shape within a few times 2^-1022 and a tilt 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).
• A different sequence: every subnormal shape, below 2.2e-308. A draw is 0 with probability 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.
• A different sequence: RTDR wherever it now works in the offset from the mode, above the tilts named above and, for 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).
• A different sequence also where the default method now takes RTDR and 0.1.1 took Algorithm 1 or 3 (see Performance): 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.

Bug fixes --- moments and mode

• 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.

Bug fixes --- density and normalising constant

• 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(50