--- title: "Theoretical background of the mhn package" output: rmarkdown::html_vignette: number_sections: true vignette: > %\VignetteIndexEntry{Theoretical background of the mhn package} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 4, dpi = 72 ) set.seed(1) ``` ```{r library, include = FALSE} library(mhn) ``` ## Scope This vignette is the theoretical companion to *Introduction to the mhn Package*. The *Introduction* shows what the four exported functions (`dmhn` / `pmhn` / `qmhn` / `rmhn`) and their moment helpers do; this note explains the mathematics and the algorithms behind them. We summarise the Modified Half-Normal (MHN) family, the Fox–Wright $\Psi$ function that appears in its normalising constant, the two sampling schemes the package draws on — Sun, Kong & Pal (2023) and Gao & Wang (2025) — and how `rmhn(method = "auto")` chooses between them. The main text keeps to what is needed to understand the results and to choose a method; the proofs are in the original papers (see References). Two appendices serve readers who want to check the implementation against those papers. Appendix A collects the numerical details: how the package keeps its results accurate where the textbook formulas lose digits, its handling of extreme parameters, and the timings behind `method = "auto"`. Appendix B lists the misprints in the papers that bear on what the package computes, and the forms it uses instead. The help page of each function states what it returns and, where it has been measured, how accurate it is. ## The MHN($\alpha$, $\beta$, $\gamma$) family The MHN distribution has support $(0, \infty)$ and density $$f(x \mid \alpha, \beta, \gamma) = \frac{2 \beta^{\alpha/2} \, x^{\alpha-1} \, \exp(-\beta x^2 + \gamma x)} {\Psi[\alpha/2,\, \gamma/\sqrt{\beta}]}, \qquad x > 0,$$ with $\alpha > 0$, $\beta > 0$, $\gamma \in \mathbb{R}$. The three parameters control three independent factors of the density kernel: | Parameter | Factor it controls | Role | |-------------|------------------------------|---------------| | $\alpha$ | $x^{\alpha-1}$ | shape | | $\beta$ | $\exp(-\beta x^2)$ | Gaussian tail | | $\gamma$ | $\exp(\gamma x)$ | exponential tilt | Three familiar distributions sit inside the family as exact special cases (Sun et al. 2023, Lemma 6): | Constraint | Reduction | |------|----------------------| | $\gamma = 0$ | $X^2 \sim \mathrm{Gamma}(\alpha/2,\, \beta)$ with rate $\beta$, i.e.\ sqrt-Gamma | | $\alpha = 1$ | Truncated normal $\mathrm{TN}(\mu = \gamma/(2\beta),\, \sigma = 1/\sqrt{2\beta},\, 0,\, \infty)$ | | $\alpha = 1,\, \gamma = 0$ | Half-normal $\lvert Z \rvert$ with $Z \sim N(0,\, \sigma^2 = 1/(2\beta))$ | `rmhn()` draws the truncated normal of the second row by the method of Robert (1995). There is also a degenerate *boundary* case, not listed in the table: Sun et al. (2023, remark after Lemma 6) note that with $\beta = 0$ and $\gamma < 0$ the law is $\mathrm{Gamma}(\alpha,\, -\gamma)$, and the MHN density tends to it as $\beta \to 0^+$. The package does not handle this limit specially because $\beta > 0$ is required by the MHN parameter space; for a Gamma draw with $\beta = 0$, use `stats::rgamma` directly. For each of the three *finite* special cases in the table above — $\gamma = 0$, $\alpha = 1$, and both together — the exported functions work through the corresponding standard law — the Gamma law of $\beta X^2$ for $\gamma = 0$, the normal law for $\alpha = 1$ — bypassing the general series-and-rejection machinery. `dmhn`, `pmhn` and `qmhn` do so at those values exactly: a shape or a tilt merely close to them is a different law and takes the general route. `rmhn` draws from the standard laws also where $|\alpha - 1|$ or the scale-free tilt $|\gamma|/\sqrt{\beta}$ is below `sqrt(.Machine$double.eps)`, which changes the law drawn by less than $8.6 \times 10^{-9}$ in total variation. Appendix A.1 gives the reasons, and the exception at very large shapes. ```{r density-shape, echo = -c(2, 8), fig.cap = "MHN density as alpha crosses the log-concavity threshold; beta = 1, gamma = 1 fixed."} x <- seq(0.001, 4, length.out = 400) op <- par(mar = c(4, 4, 2, 1)) plot(x, dmhn(x, alpha = 2, beta = 1, gamma = 1), type = "l", lwd = 2, ylim = c(0, 1.4), ylab = "f(x)", main = "MHN density at beta = 1, gamma = 1") lines(x, dmhn(x, alpha = 5, beta = 1, gamma = 1), lwd = 2, col = "steelblue") lines(x, dmhn(x, alpha = 1, beta = 1, gamma = 1), lwd = 2, col = "tomato") lines(x, dmhn(x, alpha = 0.5, beta = 1, gamma = 1), lwd = 2, col = "darkorange") legend("topright", bty = "n", lwd = 2, col = c("darkorange", "tomato", "black", "steelblue"), legend = c("alpha = 0.5 (no interior mode)", "alpha = 1 (truncated normal)", "alpha = 2 (log-concave)", "alpha = 5 (log-concave, shifted right)")) par(op) ``` The figure illustrates the qualitative claims developed in the next section. The density is log-concave whenever $\alpha \geq 1$ (Sun et al. 2023, Lemma 3a). For $\alpha > 1$ it has a finite interior mode in closed form (Sun et al. 2023, Lemma 3b); the borderline case $\alpha = 1$ reduces to the truncated normal, whose mode is $\max(0,\, \gamma/(2\beta))$ — interior when $\gamma > 0$, on the boundary $x = 0$ otherwise. For $\alpha < 1$ the density is typically monotone decreasing with a boundary singularity at $x \to 0^+$. The one exception is $\gamma > 0$ together with $\alpha > 1 - \gamma^2/(8\beta)$, where the density is not unimodal: it diverges at $x \to 0^+$, dips to a local minimum, rises to a single interior local maximum, and then decays (Sun et al. 2023, Lemma 3c, which states the condition with $\geq$; at equality the two stationary points merge into an inflection point and the density is still decreasing). At $(\beta, \gamma) = (1, 1)$ this threshold is $1 - 1/8 = 0.875$, so the orange ($\alpha = 0.5$) curve in the figure sits in the strict monotone-decreasing regime, with no interior local maximum. ## Shape and moments Three facts from Sun et al. (2023) drive everything the package does with shape, moments, and the mode. **Log-concavity is a function of $\alpha$ alone.** For $\alpha \geq 1$ the density is log-concave irrespective of $\beta$ and $\gamma$ (Sun et al. 2023, Lemma 3a); for $\alpha < 1$ it is not, because $(\log f)'' = (1 - \alpha)/x^2 - 2\beta$ is positive near $0$. This threshold is what splits the Gao & Wang (2025) RTDR sampler into its log-axis and original-axis branches in §6 below. Gao & Wang put $\alpha = 1$ itself on the original axis; the package moves it to the log axis (§6). **The mode has a closed form for $\alpha > 1$** (Sun et al. 2023, Lemma 3b): $$X_{\mathrm{mode}} = \frac{\gamma + \sqrt{\gamma^2 + 8\beta(\alpha - 1)}}{4\beta}.$$ For $\alpha = 1$ the mode is $\max(0,\, \gamma/(2\beta))$ (the truncated-normal mode, Sun et al. 2023, Lemma 6b). For $\gamma \leq 0$ this is the boundary point $0$, from which the bounded density decreases; for $\gamma < 0$ Sun et al. (2023, Lemma 3d) say that no mode exists (the statement leaves out $\gamma = 0$, which the proof in their Supplementary Section 2.7 includes), and `mhn_mode()` returns that point. For $\alpha < 1$ the density is monotone decreasing unless $\gamma > 0$ and $\alpha > 1 - \gamma^2/(8\beta)$, in which case it has a local minimum followed by an interior local maximum (Sun et al. 2023, Lemmas 3c and 3d). Lemma 3d leaves out $\gamma = 0$, where the derivative of $\log f$, $-(2\beta x^2 + 1 - \alpha)/x$, is negative, and Lemma 3c includes the boundary $\alpha = 1 - \gamma^2/(8\beta)$, where there is no local maximum (§2). The helper `mhn_mode()` implements this case split. It returns `NA` for $\alpha < 1$ without an interior local maximum, and the boundary mode 0 at $\alpha = 1$ with $\gamma \leq 0$. **Moments satisfy a two-term recurrence** (Sun et al. 2023, Lemma 2b): $$E(X^{k+2}) = \frac{\alpha + k}{2 \beta}\, E(X^{k}) + \frac{\gamma}{2 \beta}\, E(X^{k+1}).$$ Together with the closed-form ratio for $E(X)$ in terms of $\Psi$ (Sun et al. 2023, Lemma 2a), this gives every higher moment in exact arithmetic. In double precision it does not: central moments formed from raw ones cancel as soon as the standard deviation is small next to the mean, and the ratio of two $\Psi$ values loses digits under a strong positive tilt. So `mhn_mean()` and `mhn_var()` use the closed forms only where they are accurate, and elsewhere, like `mhn_skewness()` and `mhn_kurtosis()` throughout, integrate the moments numerically against the density kernel, which needs no $\Psi$ at all (Appendix A.2). The *Introduction* vignette shows the helpers in action. A useful sanity bound (Sun et al. 2023, Lemma 4c) is $\mathrm{Var}(X) \leq 1/(2\beta)$ for every $\alpha \geq 1$, $\gamma \in \mathbb{R}$. ## The Fox–Wright $\Psi$ normalising constant The normalising constant in the MHN density is $$\Psi\!\left[\frac{\alpha}{2},\, z\right] \;=\; \sum_{n=0}^{\infty} \frac{\Gamma\!\big(\tfrac{\alpha + n}{2}\big)}{n!} \, z^{n}, \qquad z = \frac{\gamma}{\sqrt{\beta}}.$$ This series is absolutely convergent for every $\alpha > 0$ and every $z \in \mathbb{R}$, and is strictly positive. It depends on $(\alpha, \beta, \gamma)$ only through the two-dimensional combination $(\alpha, z)$. A useful organisational point: **$\Psi$ matters only for some of the package's exported functions**. | Function group | Needs $\Psi$? | Why | |-------------------------------------------|:-----------------:|-----| | `dmhn`, `pmhn`, `qmhn` | yes | $\Psi$ sits in the density's denominator and in the CDF series. | | `mhn_mean`, `mhn_var` | where accurate | $E(X)$ is the ratio $\Psi[(\alpha+1)/2,\, z] / (\sqrt{\beta}\,\Psi[\alpha/2,\, z])$ and $\mathrm{Var}(X)$ follows from it, wherever that is accurate; elsewhere both come from quadrature of the kernel. | | `mhn_skewness`, `mhn_kurtosis` | no | The central moments are integrated against the kernel, whose normalising constant cancels from their ratios. | | `rmhn` (Sun or RTDR path) | no | $\Psi$ enters both proposal scale and target as a common factor that cancels. | The practical consequence is that any numerical issues in evaluating $\Psi$ — accumulated rounding in `lgamma`, or truncation error in the series — can only affect `d/p/q` and the mean and variance; they cannot leak into the sampler, apart from the weight $1/(1 + w)$ of the origin that it draws for a subnormal shape (see `?rmhn`). For evaluation, the package uses the closed forms and the series of Sun et al. (2023, Supplementary Lemmas 9 and 10) where they are accurate, and numerical integration elsewhere: - At $z = 0$ exactly $\Psi[\alpha/2, 0] = \Gamma(\alpha/2)$. At $\alpha = 1$ exactly, and at $\alpha = 2$ exactly with $\gamma > 0$, $\Psi$ has a closed form involving $\Phi$ (Sun et al. 2023, Supplementary Lemma 9(c)). A shape or a tilt merely close to these values takes the general route below. - For $0 < z \leq 20$ with $\alpha \leq 100$ the package uses the series of Sun et al. (2023, Supplementary Lemma 10). All its terms are positive, so there is no catastrophic cancellation. - Everywhere else — for $\gamma < 0$, for $z > 20$ and for $\alpha > 100$ — it integrates $\Psi[\alpha/2, z] = 2\int_0^\infty u^{\alpha-1} e^{-u^2 + z u}\,du$, the normalisation of the density in Sun et al. (2023, Lemma 1a). Past $z = 20$ the series needs hundreds of terms or more, and the quadrature is both faster and more accurate. Appendix A.3 describes the integration, and Appendix A.4 how the density, distribution and quantile functions avoid forming terms of the size of $\alpha \log x$ or $z^2/4$ under a strong tilt or at a large shape. A few results of Sun et al. (2023) that the package uses are misprinted; Appendix B lists them, with the forms the package uses. ## Sampling: Sun et al. (2023) Sun et al. (2023) propose three Accept–Reject schemes, indexed by the sign of $\gamma$ and a coarse split on $\alpha$. Two of them — Algorithm 1 and Algorithm 3 — are implemented in the package and reachable via `rmhn(method = "sun")`. Algorithm 2 is not, and this section explains the reason. ### Algorithm 1 — $\gamma > 0$, $\alpha > 1$ For each parameter triple, Algorithm 1 constructs **two** proposal distributions — a Normal with mean $\mu_{\mathrm{opt}}$ and variance $1/(2\beta)$, and a $\sqrt{\mathrm{Gamma}(\alpha/2,\, \delta_{\mathrm{opt}})}$ — each with a multiplicative envelope constant $K_{1}$ and $K_{2}$. Setup chooses whichever envelope is tighter ($K_{1} \leq K_{2}$ or the reverse). The two optima $\mu_{\mathrm{opt}}$ and $\delta_{\mathrm{opt}}$ have closed forms (Sun et al. 2023, Theorem 1), so setup is cheap. The acceptance probability is bounded below by 0.8 for $\alpha \geq 4$ (Sun et al. 2023, Theorem 2e). For $1 \leq \alpha < 4$ no uniform lower bound is proved, but Figure 2 of Sun et al. (2023), computed numerically, shows high acceptance throughout that range. At extreme shapes the package departs from the choice between the two proposals (Appendix A.5). ### Algorithm 3 — $\gamma \leq 0$ When the tilt is non-positive, Theorem 3 of Sun et al. (2023) supplies a proposal kernel based on an AM–GM-type inequality. (Sun et al. print Algorithm 3 for $\gamma < 0$, Theorem 4 and Figure 4 for $\gamma \leq 0$; at $\gamma = 0$ the package uses the square-root-Gamma closed form.) With a tuning point $m > 0$ and the exponent $r = (\beta m + |\gamma|)/(2\beta m + |\gamma|)$, the candidate is $$T \sim \mathrm{Gamma}\!\big(\alpha \cdot r,\, m(\beta m + |\gamma|)\big), \qquad X = m \cdot T^{r}.$$ The construction works for every $\alpha > 0$ (the package uses it up to $\alpha = 2^{59}$; Appendix A.5). For $\alpha > 1$, Sun et al. (2023, Theorem 4c) prove an acceptance probability of at least $1/\sqrt{2} \approx 0.707$ (the theorem is stated for $\alpha > 1$; its proof in the Supplementary Material treats $\alpha \geq 3$) at the optimal matching point and at the explicit point $m_{\star} = (\gamma + \sqrt{\gamma^{2} + 8\alpha\beta})/(4\beta)$; for $\alpha \leq 1$ the procedure still samples correctly, but no bound is proved. The matching point is initialised from a heuristic adapted from Sun et al. (2023, Section 4.2): $m_{\mathrm{init}} = \alpha^{2}/\{(1 + \alpha)\sqrt{\beta}\}$ for $\alpha \leq 1.1$ (Sun et al. print $\alpha^{2}/(1+\alpha)$, which is the same point at $\beta = 1$; the package divides by $\sqrt{\beta}$ so that the matching point scales with $X$, since $\sqrt{\beta}X \sim \mathrm{MHN}(\alpha, 1, \gamma/\sqrt{\beta})$), and a combination of the mode and the rightmost inflection point for $\alpha > 1.1$ (the inflection point is a quartic root, which the package locates by bisection). The package then applies a single Newton-Raphson refinement step in both regimes, producing $m_{\mathrm{recommend}}$, whose acceptance probability is within 0.002 of the optimum across the parameter range tabulated in Sun et al. (2023, Table 1). For this matching point the $1/\sqrt{2}$ level is observed, not proved. ### Why Algorithm 2 is not implemented Algorithm 2, stated for every $\alpha > 0$ with $\gamma > 0$ and meant for $\alpha < 1$ (Sun et al. 2023, Section 4.1.2), samples via the mixture representation $$f_{\mathrm{MHN}}(x \mid \alpha, \beta, \gamma) \;=\; \sum_{i=0}^{\infty} p_{i}\, f_{i}(x \mid \alpha, \beta), \qquad f_{i} \;\sim\; \sqrt{\mathrm{Gamma}\!\big((\alpha + i)/2,\, \beta\big)}.$$ The mixture-index sampler keeps the first $M+1$ weights and bounds the rest by a geometric envelope, which needs $M \geq M^{\dagger} = \max\{\lfloor\alpha\rfloor, \lfloor\gamma^{2}/(\varepsilon_{1}^{2}\beta)\rfloor\}$, the point beyond which successive weights shrink by at least a factor $\varepsilon_{1}$ (Sun et al. 2023, Lemma 8). With $\varepsilon_{1} = 1/2$, as in their example, $M^{\dagger} = \lfloor 4\gamma^{2}/\beta\rfloor$ for $\alpha < 1$, so for large $\gamma^{2}/\beta$ the whole scheme degrades catastrophically — Gao & Wang (2025, Table 1) report a slowdown of up to $\sim 3 \times 10^{4}\times$ against RTDR in the most extreme case $(\alpha, \gamma) = (0.01, 10000)$. The package therefore declines to implement Algorithm 2 and routes $(\alpha < 1, \gamma > 0)$ to RTDR, whose acceptance Gao & Wang (2025) bound below by $1/e$ (next section). ## Sampling: Gao & Wang (2025) RTDR Gao & Wang (2025) introduce a *Relaxed* Transformed Density Rejection (RTDR) sampler that comes with a uniform lower bound on the acceptance probability across the *entire* parameter space. Gao & Wang (2025) write the density as proportional to $x^{\lambda-1}\exp(-\alpha x^{2}-\beta x)$: their $(\lambda, \alpha, \beta)$ are the package's $(\alpha, \beta, -\gamma)$, so their $\beta < 2(\sqrt{1-2\lambda}-1)$, stated for their density scaled to their $\alpha = 1$ (their Eqs. (2) and (3)), reads $\gamma/\sqrt{\beta} > 2(1-\sqrt{1-2\alpha})$ here. Below, the symbols are the package's. Three ideas drive the construction. ### $T_{c}$-concavity Transformed density rejection (TDR) builds an envelope from tangents of $T[f]$, a monotone transform of the density, and needs optimal contact points for its tangents. Gao & Wang (2025, Section 2) review it, with the earlier literature, for the family of transforms of their Eq. (6), $$T_{c}[x] = \mathrm{sign}(c) \cdot x^{c} \quad (c > -1,\, c \neq 0), \qquad T_{0}[x] = \log x,$$ and relax its condition on the contact points (below). A density $f$ is *$T_{c}$-concave* if $T_{c}[f]$ is concave; logarithmic concavity is the $c = 0$ special case. In regions (b) and (c) below, which transform the package's RTDR implementation uses on the log-axis density turns on the sign of the tilt: for $\gamma > 0$ it takes $c = -1/2$, because $\log g$ is then convex on $(-\infty, \log(\gamma/\sqrt{\beta}/4))$, on the axis $y = \log(\sqrt{\beta}\,x)$, and a log-tangent hat would fail to dominate there, while $T_{-1/2}[x] = -x^{-1/2}$ still does. For $\gamma \leq 0$ the log-axis density is log-concave and the ordinary $c = 0$ tangent hat is used. In region (d) the sufficient condition of their Theorem 3.2 no longer holds, and the envelope of Gao & Wang (2025, Section 4) is used instead. ### Working on the log axis when $\alpha < 1$ When $\alpha < 1$ the density $f(x \mid \alpha, \gamma)$ on $(0, \infty)$ blows up at the left endpoint and is not amenable to direct envelope construction. The change of variable $y = \log x$ yields $$g(y \mid \alpha, \gamma) \;=\; \exp\!\big( \alpha y - e^{2y} + \gamma e^{y} \big), \qquad y \in \mathbb{R},$$ which is bounded and unimodal (Gao & Wang 2025, Lemma 3.3). Gao & Wang's RTDR operates on $g$ for $\alpha < 1$ and on $f$ directly for $\alpha \geq 1$; the package moves $\alpha = 1$ itself to $g$, so that it works on $g$ for $\alpha \leq 1$ and on $f$ for $\alpha > 1$ (see region (b) below). ### The four regions After scale normalisation $\sqrt{\beta} X \to X$ (so we may assume $\beta = 1$), the RTDR envelope construction splits into four regions based on $(\alpha, \gamma)$. The $\gamma$ in the table below is the normalised tilt $\gamma_{\mathrm{norm}} = \gamma/\sqrt{\beta}$; for a general $\beta$ the package's `rmhn(method = "rtdr")` rescales internally and compares $\gamma_{\mathrm{norm}}$ against the boundary $2(1 - \sqrt{1 - 2\alpha})$. | Region | $\alpha$ | $\gamma$ | Envelope target | Reference | |:------:|-----------------|----------------------|------------------------------|----------------------------------------------------| | (a) | $> 1$ $^{\dagger}$ | any | $f$ log-concave | Gao & Wang (2025, Theorem 3.1) | | (b) | $[1/2,\, 1]$ $^{\dagger}$ | any | $g$ log-concave if $\gamma \leq 0$, else $T_{-1/2}$-concave | Gao & Wang (2025, Corollary 3.1; Theorem 3.2) | | (c) | $(0,\, 1/2)$ | $\leq 2(1 - \sqrt{1 - 2\alpha})$ | $g$ log-concave if $\gamma \leq 0$, else $T_{-1/2}$-concave | Gao & Wang (2025, Lemma 3.5; Theorem 3.2) | | (d) | $(0,\, 1/2)$ | $> 2(1 - \sqrt{1 - 2\alpha})$ | $g$, with $\log g$ inflected at $y^{*} = \log(\gamma/4)$ | Gao & Wang (2025, Theorem 4.4) | $^{\dagger}$ The boundary between (a) and (b) is the package's: Gao & Wang (2025) put $\alpha = 1$ in region (a). Region (a) is the easiest: $f$ is log-concave, and two tangent lines and a constant segment suffice (Gao & Wang 2025, Eq. (11)); for $\alpha \le e$ and $\gamma/\sqrt{\beta} \le \sqrt{8\{\sqrt{(e-1)(\alpha-1)}-(\alpha-1)\}}$ the left tangent is dropped and the constant runs from 0 (their Eq. (10)). Gao & Wang put the boundary $\alpha = 1$ in region (a) (their Theorem 3.1, and `if(lambda >= 1)` in the R code of their Appendix B). The package handles it in region (b) instead, which their Theorem 3.2 also covers ($g$ is $T_{-1/2}$-concave for $\alpha \geq 1/2$): $f$ is log-concave at $\alpha = 1$ too, but the envelope of region (a) is built around an interior mode in $x$, and at $\alpha = 1$ with $\gamma \leq 0$ the density's maximum is at the boundary $x = 0$, not at an interior mode. On the log axis $g$ has a finite mode for every $\alpha > 0$ (its $x = e^{y}$ is positive), so the region (b) construction applies without a special case. Regions (b) and (c) share one construction on $g$; the sign of $\gamma$ picks the transform. Since $\mathrm{d}^{2} \log g / \mathrm{d} y^{2} = \gamma e^{y} - 4 e^{2y}$ is negative for every $y$ exactly when $\gamma \leq 0$, a $T_0$ log-tangent hat serves on that side, while for $\gamma > 0$ only $T_{-1/2}$-concavity survives and the hat switches to inverse-square pieces. Region (d) is the hardest — $\log g$ has an inflection point, and the envelope adds a chain of secant segments on the log-convex side. In the envelope of Gao & Wang (2025, Eqs. (14)-(15)), which their Theorem 4.4 validates, $K = \lceil\rho\rceil$ breakpoints sit on that side, and $\rho$ grows roughly linearly in $\gamma$ ($K = 127909$ at $\gamma = 10^{5}$ with the package's contact points; $\rho$ depends on where $t_l$ falls in its band). When the envelope has a left tangent line, the package stops the chain as soon as a breakpoint falls on or below it, the test used in the R code of their Appendix B; without one, which happens only for $\gamma$ below about 6.3, $\lceil\rho\rceil$ is at most eight. So only a handful of breakpoints is ever built, and just one once $\gamma$ exceeds about 32. The left tangent covers the breakpoints the package leaves out; their Eq. (15) misprints its interval (Appendix B.2). ### Why "relaxed" The classical TDR construction asks for *optimal* tangent contact points, which generally have no closed form and must be found numerically, by bisection or Newton–Raphson (Gao & Wang 2025, Section 2). Gao & Wang (2025, Corollary 2.1 and Eq. (7)) show that the rejection constant stays $\leq e \approx 2.718$ (equivalently, acceptance probability $\geq 1/e \approx 0.368$) as long as the contact points lie *anywhere* inside a fixed interval — namely $$\log\!\frac{f(m)}{f(t)} \in \begin{cases} [0.46,\, 2.49] & f \text{ log-concave (region (a))}, \\ [0.46,\, 2.49] & g \text{ log-concave (regions (b), (c), } \gamma \leq 0), \\ [0.93,\, 1.99] & g \text{ $T_{-1/2}$-concave (regions (b), (c), } \gamma > 0). \end{cases}$$ Setup is therefore a handful of cheap iterations that terminate as soon as the bracket is hit, with no need for a precise Newton solution. This is why Gao & Wang (2025) call the setup fast: compared with classical TDR, which must solve for optimal contact points, or with generic methods such as adaptive rejection sampling and numerical inversion, it suits varying-parameter settings such as a Gibbs sampler, where $(\alpha, \beta, \gamma)$ changes every iteration and each setup serves a single draw. They make no such claim against Sun et al. (2023), whose generator they find comparable for $\gamma < 0$ or $\alpha > 1$ (their Section 5). Algorithm 1's setup costs less than RTDR's, so `method = "auto"` keeps it for single draws with $\gamma > 0$, $\alpha > 1$, and §7 explains how it splits both signs of $\gamma$ between the Sun algorithms and RTDR by batch size and, for $\gamma > 0$, by tilt. RTDR's own advantages are a lighter per-proposal cost, which wins every large $\gamma < 0$ batch and the large $\gamma > 0$ batches at mild tilt (§7), and the $1/e$ bound below; it is also the package's only sampler for $\alpha < 1$, $\gamma > 0$. Composing the bounds of Gao & Wang (2025) — Theorem 3.1 for region (a), Corollary 2.1 with Eq. (7) for regions (b) and (c), whose concavity their Theorem 3.2 establishes, and Theorem 4.4 for region (d) — yields the headline guarantee of RTDR: for **every** $(\alpha, \beta, \gamma)$ with $\alpha > 0$, $\beta > 0$, $\gamma \in \mathbb{R}$, the acceptance probability of their construction is at least $1/e \approx 0.368$, provided its contact points lie in the band of their Eq. (7). The package's search for those points does not always reach that band (see `?rmhn`), and the acceptance measured stays above $1/e$ there too. Sun et al.'s (2023) bounds are uniform only on parts of the parameter space: at least 0.8 for Algorithm 1 when $\alpha \geq 4$ (Theorem 2e), at least $1/\sqrt{2}$ for Algorithm 3 when $\alpha > 1$ at the optimal matching point $m_{\mathrm{opt}}$ and at $m_{\star}$ (Theorem 4c; the package's one-step Newton $m$ is not covered), and no uniform one for Algorithm 2, whose bound on its index step needs $M$ to grow with $\gamma^{2}/\beta$ (§5.3). A quick sanity check that the RTDR sampler converges to the right distribution in all four regions — empirical means from RTDR vs theoretical means across regions (a), (b), (c), (d), each based on 5000 draws: ```{r rtdr-regions} cases <- data.frame( region = c("(a)", "(b)", "(c)", "(d)"), alpha = c(2.0, 0.7, 0.3, 0.3), gamma = c(0.5, 1.0, 0.5, 5.0) ) set.seed(42) empirical <- mapply( function(a, g) mean(rmhn(5000, alpha = a, beta = 1, gamma = g, method = "rtdr")), cases$alpha, cases$gamma ) theoretical <- mapply(mhn_mean, cases$alpha, beta = 1, cases$gamma) data.frame(cases, empirical_mean = round(empirical, 4), theoretical_mean = round(theoretical, 4)) ``` The four rows hit the right means to within about two Monte-Carlo standard errors — at $n = 5000$ those errors run from $0.005$ to $0.011$, and the gaps are all below $0.012$. ## The `rmhn(method = "auto")` dispatch `rmhn` exposes three values of `method`: `"sun"` (the Sun Algorithm 1 / 3 path), `"rtdr"` (the Gao & Wang RTDR path), and `"auto"` (default). The first two select the sampler used for the general case; the closed-form special cases below are taken first whatever `method` is set to, since they are exact and need next to no setup. The `"auto"` route is the package's recommendation; it combines the special cases and the three samplers into one decision tree, whose ordinary thresholds come from timing the samplers against each other: ``` rmhn(n, alpha, beta, gamma, method = "auto") │ ▼ ┌─ alpha = 1, gamma = 0 → half-normal (Sun et al. 2023, Lemma 6c) │ ├─ gamma = 0 → sqrt-Gamma (Sun et al. 2023, Lemma 6a) │ (above alpha = 2^60: the Normal proposal of Sun Algorithm 1) │ ├─ alpha = 1 → truncated normal (Sun et al. 2023, Lemma 6b) │ ├─ alpha < 1 and gamma > 0 │ → RTDR (regions (b) / (c) / (d), Gao & Wang 2025, Theorems 3.2, 4.4) │ ├─ alpha > 1 and gamma > 0 │ ├─ alpha ≤ 2^59, gamma / sqrt(beta) ≤ 2 sqrt(alpha), │ │ samples_per_setup ≥ N+ → RTDR (region (a)) │ └─ otherwise → Sun Algorithm 1 │ (Sun et al. 2023, Theorems 1, 2) │ where N+, for gamma / sqrt(beta) ≤ sqrt(alpha) / above it, is │ 5 / 25 (alpha < 10), 25 / 50 (10 ≤ alpha < 1e3), │ 50 / 100 (1e3 ≤ alpha < 1e8), 100 / 1000 (alpha ≥ 1e8), │ and 5 / 25 where RTDR's setup starts next to the mode (*) │ └─ gamma < 0 ├─ alpha > 2^59 → RTDR ├─ -gamma / sqrt(beta) > 1e150 → Sun Algorithm 3 ├─ alpha < 2.5e-324 * (-gamma / sqrt(beta)) │ → Sun Algorithm 3 (the log-axis mode rounds to 0) ├─ samples_per_setup ≥ N*, or RTDR's setup │ starts next to the mode (*) → RTDR (regions (a)-(c)) └─ otherwise → Sun Algorithm 3 where N* = 25 for alpha < 0.1, 5 for 0.1 ≤ alpha < 1, 1 for 1 < alpha < 1e3, 10 for 1e3 ≤ alpha < 1e5, 25 for 1e5 ≤ alpha < 1e8, 50 for alpha ≥ 1e8 (*) region (a), alpha > 1: gamma / sqrt(beta) < -1e10 sqrt(alpha), or the centred construction, from a shape of about 2.4e12 to 3.1e12 at moderate tilts, of about 2.3e13 near a tilt of -alpha / e and of about 4.5e13 at any tilt; the default method asks the setup's own test, and these shapes are approximate alpha = 1 and gamma = 0 above are taken to within sqrt(.Machine$double.eps), as described in Section 2 ``` Write $S$ for the `samples_per_setup` quantity that the branches above compare with $N^{+}$ and $N^{*}$, $$S \;=\; \begin{cases} n, & L_{\alpha} = L_{\beta} = L_{\gamma} = 1, \\ \lfloor n / P \rfloor, & 2 \le P \le 256 \text{ and } 4P \le n, \\ 1, & \text{otherwise}, \end{cases}$$ where $L_{\alpha}$, $L_{\beta}$, $L_{\gamma}$ are the lengths of the three parameter vectors and $P = \operatorname{lcm}(L_{\alpha}, L_{\beta}, L_{\gamma})$ is the period after which the recycled triples repeat. In the middle case `rmhn` keeps a setup for each position within the period, and each serves $\lfloor n/P \rfloor$ draws; otherwise consecutive draws may need different setups, and a call with parameter vectors of length $n$ counts one draw per setup, as the loop of single draws in a Gibbs sampler does. The intuition follows from decomposing the expected time per draw as $$E[\text{time}/\text{draw}] \;=\; \frac{T_{\text{setup}}}{S} \;+\; C_{\text{reject}} \times T_{\text{per-proposal}},$$ where $T_{\text{setup}}$ is the one-off cost of preparing the proposal (the $m_{\text{init}}$ heuristic plus at most one Newton step for Sun Algorithm 3; an iterative contact-point search for RTDR), $C_{\text{reject}} = 1/\text{acceptance}$ is the expected number of proposals per accepted draw, and $T_{\text{per-proposal}}$ is the cost of one proposal cycle (sample plus acceptance log-ratio evaluation). For $\gamma < 0$ and $\alpha < 1$ the two samplers sit on opposite ends of this trade-off: - **Sun Algorithm 3** has a *cheap setup* — $m_{\text{init}}$ is explicit for these shapes, followed by a single Newton refinement — but a *more expensive per-proposal*: every proposal calls `rgamma` and evaluates the fractional power $X = m \cdot T^{r}$. - **RTDR** (regions (b) and (c) for these shapes) has a *more expensive setup*: each contact point $t_{l}$, $t_{r}$ is found by an iterative search that terminates as soon as the log-ratio of the density at the mode to the density at the contact point lands in the bracket $[0.46,\, 2.49]$. But the *per-proposal* is light: a draw from the piecewise envelope and the acceptance test need only three uniform variates and a few `log` and `exp` evaluations. When $S$ is large the per-proposal term dominates and RTDR wins; when it is small — the Gibbs case where each iteration redraws $(\alpha, \beta, \gamma)$ — the setup term does. The crossovers in the tree come from timing the samplers against each other on an Apple M1 Ultra, and they move from machine to machine; Appendix A.6 gives the timings and the margins. The branches that start at extreme values — shapes above $2^{59}$, tilts beyond $10^{150}$, and laws RTDR cannot place — come from the limits of double precision instead, and ordinary parameters never reach them; Appendix A.5 explains them. A small concrete demonstration: with `set.seed` aligned, `method = "auto"` and the *expected* forced method produce bit-identical draws on each branch. ```{r auto-dispatch} # Branch: alpha < 1 and gamma > 0 --> expects RTDR set.seed(7); auto1 <- rmhn(20, alpha = 0.7, beta = 1, gamma = 2) set.seed(7); rtdr1 <- rmhn(20, alpha = 0.7, beta = 1, gamma = 2, method = "rtdr") identical(auto1, rtdr1) # Branch: alpha > 1, gamma > 2 sqrt(alpha beta) --> expects Sun Algorithm 1 set.seed(7); auto2 <- rmhn(20, alpha = 3, beta = 1, gamma = 5) set.seed(7); sun2 <- rmhn(20, alpha = 3, beta = 1, gamma = 5, method = "sun") identical(auto2, sun2) # Branch: alpha > 1, gamma <= sqrt(alpha beta), 5 or more draws per setup # --> expects RTDR set.seed(7); auto3 <- rmhn(20, alpha = 3, beta = 1, gamma = 1.5) set.seed(7); rtdr3 <- rmhn(20, alpha = 3, beta = 1, gamma = 1.5, method = "rtdr") identical(auto3, rtdr3) # Branch: gamma < 0, 0.1 <= alpha < 1, fewer than 5 draws per setup # --> expects Sun Algorithm 3 set.seed(7); auto4 <- rmhn(4, alpha = 0.5, beta = 1, gamma = -1) set.seed(7); sun4 <- rmhn(4, alpha = 0.5, beta = 1, gamma = -1, method = "sun") identical(auto4, sun4) ``` The four `TRUE` outputs are the same invariant that the package's tests enforce as a regression test. ## References Sun, J., Kong, M., & Pal, S. (2023). The Modified-Half-Normal distribution: Properties and an efficient sampling scheme. *Communications in Statistics - Theory and Methods*, 52(5), 1591--1613. Gao, F., & Wang, H.-B. (2025). Generating modified-half-normal random variates by a relaxed transformed density rejection method. *Communications in Statistics - Simulation and Computation*, 1--25. Advance online publication. Robert, C. P. (1995). Simulation of truncated normal variables. *Statistics and Computing*, 5(2), 121--125. ## Appendix A. Numerical details {-} This appendix is for readers who want to check the implementation against the papers. It describes how the package keeps its results accurate where the textbook formulas lose digits, how it treats parameters at the edges of double precision, and the timings behind `method = "auto"`. None of it is needed to use the package; the help pages state what each function returns and, where it has been measured, how accurate it is. ### A.1 Special cases and their neighbourhoods {-} `dmhn`, `pmhn` and `qmhn` use the standard laws of §2 only at the special values exactly. A shape or a tilt merely close to them is a different law: at $\alpha = 1 + \delta$ the log density moves by about $\delta(\log(\sqrt{\beta}x) - E\log U)$, $U = \sqrt{\beta}X$, which next to $x = 10^{-300}$ is about $10^{-5}$ for $\delta = 1.5 \times 10^{-8}$, so such a triple takes the general route, which holds its accuracy up to the special values. `rmhn` draws from the standard laws also where $|\alpha - 1|$ or the scale-free tilt $|z|$, $z = \gamma/\sqrt{\beta}$, is below `sqrt(.Machine$double.eps)`. The law drawn is then within a total variation distance of about $|\delta|\,E|\log U - E\log U|/2$ or $|z|\,E|U - EU|/2$ of the one requested, and of $E|\delta L + zV|/2$, with $L$ and $V$ those two centred terms, where both are that small and the half-normal is drawn; this is below $8.6 \times 10^{-9}$. The band also spares `method = "sun"` a case it cannot draw: of Sun et al.'s samplers only Algorithm 2, which the package does not implement (§5.3), covers $\alpha < 1$ with $\gamma > 0$, a case `method = "sun"` would otherwise meet next to $\alpha = 1$ from below under a positive tilt, and next to $\gamma = 0$ from above for any $\alpha < 1$. The one exception to the use of the standard laws is a zero or negligible tilt at a very large shape. The Gamma law then has a relative width of about $\sqrt{2/\alpha}$, and the rounding of $\beta x^2$, or of R's `rgamma`, would move it by a visible part of that width. So `dmhn`, `pmhn` and `qmhn` use it only up to $\alpha = 10^6$, beyond which they treat $\gamma = 0$ like any other tilt (see Appendix A.4), and `rmhn` only up to $\alpha = 2^{60}$, beyond which a rejection sampler placed at the mode takes over (Appendix A.5). ### A.2 Moments and mode {-} Throughout, $z = \gamma/\sqrt{\beta}$. The recurrence of §3, together with the closed-form ratio for $E(X)$ in terms of $\Psi$ (Sun et al. 2023, Lemma 2a), gives every moment in exact arithmetic. In double precision it does not: central moments formed from raw ones cancel as soon as the standard deviation is small next to the mean, and the ratio of two $\Psi$ values whose logarithms are each about $z^2/4$ loses digits under a strong positive tilt. The helpers therefore use the closed forms only where they are accurate, and elsewhere integrate the moments numerically against the density kernel — in the offset from its peak where it has one, and on a logarithmic axis out from the origin where it has none — which needs no $\Psi$ at all. **Mean.** `mhn_mean()` forms the ratio $E(X) = \Psi[(\alpha+1)/2,\, z] / (\sqrt{\beta}\,\Psi[\alpha/2,\, z])$ from the logarithms of the two $\Psi$ values. It uses the ratio for $z \le 12$ as long as these two logarithms are small enough to leave it accurate. For a larger positive tilt, or a very large $\alpha$, both logarithms are large while their difference is not, and the mean is integrated instead. For $z > 12$ and a shape so small that most of the mass lies by the origin, that integration holds the mean as a small difference; where it finds the mean below a thousandth of the peak it is centred on, the ratio is used instead. Under a strong negative tilt, $z < -30$, the ratio serves only for $\alpha < 6$, and there it is capped at the upper bound $(\gamma + \sqrt{\gamma^2 + 8\alpha\beta})/(4\beta)$ of Sun et al. (2023, Lemma 4a), which holds for every $\alpha > 0$ by the remark after the lemma; from $\alpha = 6$ on, the mean is integrated there as well. At $\gamma = 0$ the ratio reduces to $\Gamma((\alpha+1)/2) / \big(\sqrt{\beta}\, \Gamma(\alpha/2)\big)$ (Sun et al. 2023, Lemma 6a). For $\alpha \ge 1$ it is computed as $\sqrt{\pi} / \big(\sqrt{\beta}\, B(\alpha/2,\, 1/2)\big)$, with $B$ the beta function, whose logarithm R's `lbeta` forms without the cancellation between the two log-gamma values that loses digits at a large $\alpha$. **Variance.** `mhn_var()` uses the form of Sun et al. (2023, Lemma 2c), $\mathrm{Var}(X) = \alpha/(2\beta) + E(X)\,\big(\gamma/(2\beta) - E(X)\big)$, where it is accurate. Its terms cancel under a strong tilt of either sign, where the form also magnifies any error in $E(X)$ by about $z^2/2$; at a large $\alpha$ it magnifies that error by about $4\alpha$. Wherever the combined error, from the rounding of its terms and from the error of the mean, would exceed about $10^{-8}$ of the result, the variance is integrated instead. This is so for $z > 12$, where the mean itself is integrated, for $z < -30$, and for $\alpha \ge 1000$ unless $\gamma = 0$. **Skewness and kurtosis.** `mhn_skewness()` and `mhn_kurtosis()` integrate the central moments. The raw moments of the Lemma 2b recurrence, started from the mean of Lemma 2a, serve only where the integration cannot be set up, as for a shape so small, under a positive tilt, that the mean falls far below the peak on which the integration is centred; the standard deviation of such a law is many times its mean, and the raw moments do not cancel there. **Subnormal shapes.** For a shape below the smallest normal double, $2.2 \times 10^{-308}$, the moments are taken from the forms in which the shape only multiplies, the same forms as the density takes there (Appendix A.4), with $w$ the mass at the peak of the kernel against the mass by the origin. While $w$ is well below 1, up to a tilt $z$ of about 53 at the largest subnormal shape and 54.5 at the smallest, the skewness grows like $1/\sqrt{\alpha}$ and the excess kurtosis like $1/\alpha$ as the shape falls. Where $w$ is about 1, for $z$ from about 53.3 to 54.7, the skewness turns negative and in the kurtosis the peak takes over. Past the tilt where the forms hand over, 55 to 56.3, neither depends on $\alpha$. **Mode.** The mode is a closed form in every case (§3), and what matters is how it is evaluated. For $\alpha > 1$ it is the positive root of $2\beta x^2 - \gamma x - c = 0$, $c = \alpha - 1$. The printed form $(\gamma + \sqrt{\gamma^2 + 8\beta c})/(4\beta)$ subtracts two nearly equal quantities when $\gamma < 0$, since $\sqrt{\gamma^2 + 8\beta c}$ tends to $|\gamma|$ as $|\gamma|$ grows; multiplied through by the conjugate it becomes $2c/(\sqrt{\gamma^2 + 8\beta c} - \gamma)$, which only adds for $\gamma < 0$, and each sign takes the form that only adds. Where the squares and products of these forms would overflow, or lose digits below the smallest normal double, while the root itself lies well inside the range of doubles, the root is formed from ratios instead: with $r = \sqrt{8\beta c}$, the larger of $|\gamma|$ and $r$ is factored out, leaving a quantity within a factor $1 + \sqrt{2}$ of the root times one of order one. The interior maximum of Sun et al. (2023, Lemma 3c), for $\alpha < 1$ and $\gamma > 0$, is treated alike. The truncated-normal mode $\max(0, \gamma/(2\beta))$ is taken at $\alpha = 1$ itself only, with no tolerance band around 1: the root of Lemma 3b tends to it as $\alpha \to 1^{+}$, and for $\gamma > 0$ so does the root of Lemma 3c from below. ### A.3 Evaluating $\Psi$ {-} For evaluation, the package uses the closed forms and the series of Sun et al. (2023, Supplementary Lemmas 9 and 10) where they are accurate, and numerical integration elsewhere: - At $z = 0$ exactly $\Psi[\alpha/2, 0] = \Gamma(\alpha/2)$. A tilt merely close to 0 takes the general route below. - At $\alpha = 1$ exactly, and at $\alpha = 2$ exactly with $\gamma > 0$, $\Psi$ has a closed form involving $\Phi$ (Sun et al. 2023, Supplementary Lemma 9(c)). For $\alpha = 1$ under a strong negative tilt ($z \leq -4$) it is written through the scaled complementary error function, so that the $z^2/4$ in its exponent cancels analytically. A shape merely close to 1 or 2 takes the general route below. - For $0 < z \leq 20$ with $\alpha \leq 100$ the package uses the series of Sun et al. (2023, Supplementary Lemma 10), truncated at an upper bound on the point that lemma gives, or earlier once the terms left can no longer change the sum. All its terms are positive, so there is no catastrophic cancellation. - Everywhere else — for $\gamma < 0$, for $z > 20$ and for $\alpha > 100$ — it integrates $\Psi[\alpha/2, z] = 2\int_0^\infty u^{\alpha-1} e^{-u^2 + z u}\,du$, the normalisation of the density in Sun et al. (2023, Lemma 1a). Past $z = 20$ the series needs hundreds of terms or more, and the quadrature is both faster and more accurate. Above $\alpha = 100$ the series is no more accurate than the quadrature and, except at the smallest tilts, where the two take about the same time, slower. The integral is taken in one of two variables. In either, the integrand is divided by its value at its peak, split there, and cut off where it has fallen by $L = -\log(\mathrm{tol}) + 40$ below the peak, about 58 at the working tolerance; the panels go to Boost.Math quadrature rules. - **In $v = \log u$**, for $z \leq 0$, for $z > 10^3$, for $\alpha > 10^6$, and for a shape below 6 whose peak lies within $\sqrt{L} \approx 7.6$ of the origin. On this axis the integrand is $u^{\alpha} e^{-u^2 + z u}$, which for either sign of $z$ has a single analytic peak, at the positive root $u^{*}$ of $2u^2 - z u - \alpha = 0$. With $u = u^{*} e^{w}$, $$\Psi[\alpha/2, z] = 2\, (u^{*})^{\alpha} e^{-(u^{*})^2 + z u^{*}} \int_{-\infty}^{\infty} \exp\!\big(-\alpha(e^{w} - 1 - w) - (u^{*})^2 (e^{w} - 1)^2\big)\,dw.$$ Both terms of the exponent are at most 0 and vanish at the peak, so nothing large is subtracted however strong the tilt, and neither the origin nor the power $u^{\alpha-1}$, whose derivatives are unbounded there, appears. The right side is one 31-point Gauss–Kronrod panel. Far out, the left side falls only like $\alpha |w|$ — that is the mass $u^{\alpha-1}$ keeps near the origin — and it is taken on panels $[-2d, -d]$, doubling outwards, until a bound on what remains falls below $e^{-L}$ of the sum, or until a head next to the origin is reached. That head, below a point $x_0$ at which $e^{-u^2 + z u}$ agrees with $1 + z u - u^2$ to within $10^{-17}$, is integrated in closed form; it holds most of the mass when $\alpha$ is tiny — below any representable abscissa once $\alpha$ is below about $10^{-3}$. - **In $u$**, under a positive tilt with $z \leq 10^3$ and $\alpha \leq 10^6$ otherwise. The integrand then has an interior peak, near $z/2$ under a strong tilt, and is a bell about one unit wide there. One 31-point Gauss–Kronrod panel of at most $\sqrt{L}$ units on each side of the peak integrates it without subdividing (for $\alpha < 1$ the right panel goes to the double-exponential `tanh_sinh` rule instead). The range further below the peak is negligible at the working tolerance and is left out, except for a peak nearer the origin than $\sqrt{L}$ at $\alpha \geq 6$, where the left panel reaches down to the same closed-form head. In $u$ the terms $\alpha \log u$, $u^2$ and $z u$ of the exponent, of size $z^2/4$ and more at the peak, round by more than a panel's tolerance from $z$ of a few thousand and from $\alpha$ near $10^7$, and the switch to $\log u$, which forms none of them, comes before that. The density, distribution and quantile functions need $\Psi$ only through the difference between its logarithm and the kernel exponent $(\alpha - 1)\log x - \beta x^2 + \gamma x$. Under a strong positive tilt both are about $z^2/4$ while that difference is of order one, so for $z > 10^3$ both are formed with $z^2/4$ removed analytically. Where even that does not suffice, at a very large shape or tilt, those functions measure the law from its peak instead, with the integral over $w$ above as its normalisation (Appendix A.4). ### A.4 Density, distribution and quantile functions {-} #### Shared by the three functions {-} Throughout, $z = \gamma/\sqrt{\beta}$ is the scale-free tilt and $m$ is the larger of $\alpha$ and $z$. The three functions evaluate the law in one of two ways. The *ordinary forms* are the logarithm of the density formula, the series of Sun et al. (2023, Lemma 1b) for the distribution function, and the quadratures of the density. They are used only while $m \le 10^6$ and $m(1 + |\log v^{*}|)$, the size of their largest terms, is at most $10^7$. At $\beta = 1$ this allows a shape up to $10^6$ or a tilt up to about $7.3 \times 10^5$; at $\beta = 10^{-300}$ or $\beta = 10^{300}$ it allows up to about $3 \times 10^4$ for either. For a shape or a tilt above $10^6$, and wherever the terms of the ordinary forms grow past that bound, the law is *measured from the peak* $v^{*}$ of the density of $\log X$, the positive root of $$2\beta x^2 - \gamma x - \alpha = 0,$$ which lies at or above the median. With $x = v^{*} e^{w}$, the kernel $x^{\alpha - 1} e^{-\beta x^2 + \gamma x}$ divided by its value at the peak, taken against $dx$, becomes $v^{*}$ times $$\exp\!\big(-\alpha(e^{w} - 1 - w) - \beta v^{*2}(e^{w} - 1)^2\big)$$ against $dw$. Both terms of the exponent are at most 0 and vanish at the peak, so nothing of the size of $\alpha \log x$, $(\alpha/2)\log\beta$ or $\gamma^2/(4\beta)$ is formed. The integral of this function over the whole line, times $v^{*}$, is the peak-normalised integral of the law, which normalises the density and both tails. $v^{*}$ is held as a double next to it plus its offset from that double, so that the offset $d = x - v^{*}$ of a point from the peak is exact, even where the law is narrower than the spacing of the doubles near its peak. The integral over $w$ is the one of Appendix A.3 with $u^{*} = \sqrt{\beta}\,v^{*}$. A subnormal shape, below the smallest normal double $2.2 \times 10^{-308}$, is computed apart, at $\gamma = 0$ too. There half of $\alpha$ is not always a double and $1/\alpha$ can overflow, while $t^{\alpha} = 1$ to double precision at every positive double $t$. The law then takes a form in which $\alpha$ only multiplies. Its density is $$f(x) = \frac{\alpha\, e^{-\beta x^2 + \gamma x}}{x\,(1 + w)}, \qquad w = \alpha K_0(z),$$ with $$K_0(z) = \int_0^1 \frac{e^{-t^2 + z t} - 1}{t}\,dt + \int_1^\infty \frac{e^{-t^2 + z t}}{t}\,dt,$$ and $1/(1 + w)$ is the probability that $X$ lies below the smallest positive double. `pmhn()`, `qmhn()`, `rmhn()` and the moment functions take the same form. For a large tilt $w$ is the mass at the peak of the kernel against the mass by the origin; it is below $10^{-200}$ in size for $z \le 30$, and negative below $z$ of about 0.3. Below a tilt of about 52 (53 at the smallest shapes) $1/(1 + w)$ is 1 in double precision, with the consequences for `pmhn()`, `qmhn()` and `rmhn()` that the help page of `dmhn()` lists. The form serves the whole law 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}$, about $5.4 \times 10^{-20}$, of the whole, and the density, the upper tail, the moments and the law of the draws no longer depend on $\alpha$; the lower tail far below the peak, of which that mass is more than $2^{-64}$, is still taken from the same form up to where it reaches one half. #### Density {-} The density is computed in log space, to avoid underflow and overflow. At $\gamma = 0$ exactly, the closed form of the square-root-Gamma law is used only for $\alpha \le 10^6$. Above that the Gamma law of $\beta X^2$ is so narrow that rounding $\beta x^2$ moves it by a visible part of its width, by up to about $10^{-5}$ of the log density at $\alpha = 10^{20}$, and $\gamma = 0$ is measured from the peak like any other tilt. At $\alpha = 1$ exactly the law is the normal law of mean $\mu = \gamma/(2\beta)$ and standard deviation $\sigma = 1/\sqrt{2\beta}$ restricted to $(0, \infty)$. Write $z_0 = -\gamma/\sqrt{2\beta}$ for the truncation point 0 in standard units. Where 0 lies within 1000 standard deviations of the mean, $|\gamma|/\sqrt{2\beta} \le 1000$, the density is `dnorm` divided by its upper tail at 0. Farther out the same ratio is written in one of two other ways. Where 0 lies above the mean, $z_0 > 10^3$, it is written with the normal hazard $\phi(z_0)/\bar\Phi(z_0)$, so that two logarithms of about $-\gamma^2/(4\beta)$ are not subtracted. Where 0 lies below the mean, $z_0 < -10^3$, it is written from the standardised point $(x - \mu)/\sigma$, formed directly: the offset $2\beta x - \gamma$ is formed in one rounding, as $\beta x - \gamma/2$, and so also where $2\beta$ overflows. In this way the density holds to the Exponential and normal limits it tends to. In the ordinary form the log density adds terms of size $\alpha \log x$ and $(\alpha/2)\log\beta$, and for $z > 10^3$ measures $x$ from the centre $\gamma/(2\beta)$ rounded to a double. Their rounding passes $10^{-7}$ at $\alpha = 10^8$, and $10^{-6}$ at $z = 10^{10}$ unless $\beta$ is a power of 2. The bound of $10^7$ on the size of the largest terms (see the shared section) keeps that rounding to about $2 \times 10^{-9}$. Measured from the peak, with $d = x - v^{*}$ and $\rho = d/v^{*}$, the log density is $$\log f(x) = (\alpha - 1)\big(\log(1 + \rho) - \rho\big) - \rho - \beta d^2 - \log\!\Big(v^{*} \int_{-\infty}^{\infty} \exp\!\big(-\alpha(e^{w} - 1 - w) - \beta v^{*2}(e^{w} - 1)^2\big)\,dw \Big),$$ the last term being the logarithm of $v^{*}$ times the peak-normalised integral of the law. Neither of the large terms of the ordinary form is formed, and $d$ is exact. The error is then about $10^{-12}$, relative to the size of the log density where that is large. Together, the two bounds give the range stated on the help page: for $\alpha \le 10^{15}$, $|z| \le 10^{15}$ and $\beta$ in $[10^{-300}, 10^{300}]$ the log density is accurate to the working tolerance `sqrt(.Machine$double.eps)`. #### Distribution function {-} The series on the help page is that of Sun et al. (2023, Lemma 1b) written in the package's notation. There each term carries the factor $2\beta^{\alpha/2}\gamma^{i}\beta^{-(\alpha+i)/2}/(2\, i!)$, which reduces to $z^{i}/i!$, and the paper's lower incomplete gamma function $\gamma(s, y)$ equals $\Gamma(s)\,P(s, y)$. The infinite sum is truncated at an upper bound on the length $K = \max\{K_1, K_2\}$ of Sun et al. (2023), Supplementary Lemma 10(d), obtained from the geometric decay of the terms that Lemma 10(a, b) establish; this makes the truncation residual bounded by the working tolerance `sqrt(.Machine$double.eps)` divided by $\Psi$. For $\gamma > 0$, where every term is positive, the sum is truncated earlier, where that bound is below the rounding of the sum, so that the residual is below the larger of that and $2^{-53} F$. The series is tried first for $\alpha < 1$ at every tilt up to $z = 1000$, and for $\alpha \ge 1$ only where $|z| \le 20$ and $\alpha \le 100$. There a numerical integration of the density on $[0, q]$ replaces it in three cases: for $\gamma < 0$, where the rounding of the alternating terms, set against the cancellation between them, could make its relative error exceed `sqrt(.Machine$double.eps)`; where it would need more than $10^6$ terms, which happens for $|z|$ above about 830, and so only for $\alpha < 1$; and where $\beta q^2$ is below the smallest normal double. Elsewhere, outside the ranges measured from the peak, the integration is used directly, being the faster of the two there. It is Gauss–Kronrod for $\alpha \ge 1$. For $\alpha < 1$, where the density is unbounded at the origin, it is a closed form next to the origin, then tanh-sinh up to $q$ or up to the interior peak that a positive tilt can give the density, whichever comes first, and Gauss–Kronrod past that peak. For $z > 1000$ the integral and $\Psi$ are both taken with $\gamma^2/(4\beta)$ removed from their logarithms, which are otherwise too large for their difference to be accurate. The upper tail $P(X > q)$ is computed by Gauss–Kronrod quadrature of the density over $(q, \infty)$ on a logarithmic axis. Each tail is computed directly on its own side of the peak $v^{*}$, and above the peak the lower tail is returned as one minus the upper tail. With `log.p = TRUE` both tails are formed on the log scale throughout. At or below the median $\log P(X > q)$ is about $-F$, and the logarithm of the quadrature carries an absolute error of about $\epsilon$ times the size of the terms it is assembled from, $$|\log\Psi| + (\alpha/2)\,|\log\beta| + m\,(1 + |\log v^{*}|),$$ with $\gamma^2/(4\beta)$ removed from $\log\Psi$ for $z > 1000$, as in the quadratures. Wherever that could pass $10^{-9}$ of $F$, the log upper tail is formed instead as $\log(1 - F)$ from $\log F$, which keeps the relative accuracy of $F$; for parameters of order one this happens only below $F$ of about $10^{-5}$. For $\alpha < 1$ with $|z| > 20$, where the series costs milliseconds, $F$ is taken for this from the numerical integration wherever that agrees with the upper tail to within the error of the latter. The lower tail has the matching switch. Below the peak $F$ can still exceed one half, and for a small shape it is close to 1 over most of the range below the peak, one less a tail proportional to $\alpha$. There $\log F$, about $-P(X > q)$, is assembled from terms of the size of $\log\Psi$, $(\alpha/2)\log\beta$ and $(\alpha - 1)\log q$ that cancel, and carries an absolute error of about $\epsilon$ times their size. Wherever that could pass $10^{-9}$ of $\log F$, and the error of the upper tail, carried into $\log F$, is estimated at less than half as much, $\log F$ is formed instead as $\log(1 - P(X > q))$ from the upper tail, which keeps the relative accuracy of $P(X > q)$. Without that switch the error would pass $10^{-9}$ of $\log F$ for a shape below about $10^{-6}$ under a small tilt, and below a smaller shape under a large one. For a subnormal shape the form of the shared section keeps the relative accuracy of $\log F$ as well. At $\gamma = 0$ exactly, $F(q) = P(\alpha/2, \beta q^2)$ is taken from `pgamma`, or, where $\beta q^2$ is below the smallest normal double, as its leading term $(\beta q^2)^{\alpha/2}/\Gamma(\alpha/2 + 1)$, which is exact to double precision there. This is done for $\alpha \le 10^6$. Above it the Gamma law is so narrow that rounding $\beta q^2$ moves it by a visible part of its width, by up to about $10^{-6}$ of $F$ at $\alpha = 10^{20}$, and $\gamma = 0$ is measured from the peak. At $\alpha = 1$ exactly the truncated-normal distribution function is taken from `pnorm`, or through the scaled complementary error function where the truncation point 0 lies three or more standard deviations above the mean $\gamma/(2\beta)$ of the untruncated normal, $z_0 = -\gamma/\sqrt{2\beta} \ge 3$. Near the truncation point the lower tail is summed from a short series, so that it keeps its relative accuracy however small it is. For $z_0$ far beyond 3 the same forms hold up to $z_0 = \infty$, where the law is Exponential with rate $|\gamma|$. For $z_0$ below $-10^3$ the standardised offset is formed from $2\beta q - \gamma$ in one rounding (as $\beta q - \gamma/2$, so also where $2\beta$ overflows), not from a rounded mean $\gamma/(2\beta)$. The lower tail within $|z_0|/2$ standard deviations of 0, where that offset is resolved only to about $\epsilon |z_0|$ of one, is taken from the ratio of the two normal lower tails written with the exact distance $q\sqrt{2\beta}$ from the truncation point. Measured from the peak, both tails are integrals in the variable $w = \log(x/v^{*})$. The tail on the side of $\log(q/v^{*})$ away from the peak is the integral of $\exp\!\big(-\alpha(e^{w} - 1 - w) - \beta v^{*2}(e^{w} - 1)^2\big)$ over that side, taken by Gauss–Kronrod quadrature in offsets from $q$, divided by the same integral over the whole line; the other tail is its complement. Nothing of the size of $\alpha \log x$ or $\gamma^2/(4\beta)$ is formed, whose rounding would otherwise leave the tails $10^{-7}$ off at $\alpha = 10^8$ and farther off as $z$ grows, and no abscissa is rounded to the doubles next to the peak, which are a fraction $z\,2^{-53}$ of the law's standard deviation apart. In the ordinary forms the largest terms are $\alpha \log x$, $(\alpha/2)\log\beta$ and the rounding of quadrature abscissae placed in $\log x$, all of size about $m(1 + |\log v^{*}|)$. The bound of $10^7$ on that size keeps their error within 10 standard deviations of the peak below $10^{-8}$; measured from the peak the error is about $10^{-14}$. Farther out in the upper tail the error of the ordinary forms grows with the distance from the peak, and close to that bound it reaches about $3 \times 10^{-8}$ at 36 standard deviations, where the tail is of order $10^{-284}$. This is the range stated on the help page: for $\alpha \le 10^{15}$, $|z| \le 10^{15}$ and $\beta$ in $[10^{-300}, 10^{300}]$ both tails are accurate to the working tolerance within 10 standard deviations of the peak, and to about twice that beyond. #### Quantile function {-} In the general case $q = F^{-1}(p)$ is found by the TOMS 748 root-finder, applied to the distribution function computed as in `pmhn()`, or, for `lower.tail = FALSE`, to the logarithm of its directly computed upper tail. The initial bracket is $[s\sqrt{\epsilon},\, s]$ with $s = E(X) + 8\sqrt{\mathrm{Var}(X)}$, which scales with the distribution. Its left end is halved (down to the smallest positive double) and its right end doubled (up to the largest double) until it brackets the target probability, each step keeping the bracket a factor of two wide. Where the root-finder does not reach its tolerance on that bracket, as for a tiny lower-tail $p$ with a large $\alpha$, the equation is solved again for $\log q$ on the logarithm of the distribution function; where that fails too, or the bracket lies near the largest double, the bracket is bisected. A lower-tail `log.p` probability below `log(.Machine$double.xmin)`, about $-708.4$, whose exponential is no longer a normal double, is solved on the logarithm of the distribution function. Near 1 the logarithm of the upper tail is about $-(1 - p)$, and its error, set by the size of the terms it is assembled from (see Distribution function above), is a relative error in $1 - p$. Where that could pass $10^{-9}$, which happens within about $3 \times 10^{-6}$ of 1 at the least, and further out for a large shape, tilt or $|\log\beta|$, the distribution function is solved for one minus the probability instead, which is exact. Likewise a lower-tail $p$ close to 1, where the distribution function is resolved only to the rounding of 1, is solved on the logarithm of the upper tail for one minus it, which is exact. That is done within $10^{-7}$ of 1 at least, and in general wherever the error of $F$ near 1 could pass $10^{-9}$ of $1 - p$ and the upper tail resolves it more accurately: $F$ is resolved there only to the rounding of 1, and, for a small shape, where it is close to 1 over most of the range, to the error of the lower tail. At parameters of order one that is so within about $10^{-5}$ of 1, and farther for a small or large shape, a large tilt or an extreme $\beta$. Where `pmhn()` takes the distribution function from the series of Sun et al. (2023), which can cost far more than the upper tail, the root of the upper tail is refined instead by Newton steps on the logarithm of the distribution function, taken as in `pmhn()` for the log upper tail. Within $10^{-7}$ of 1, where the upper tail may not resolve one minus the probability at all, that root is taken for a probability it does resolve where needed, and the steps start from it. At $\gamma = 0$ exactly, $y$ in $q = \sqrt{y/\beta}$ is taken from `qgamma`, or from the leading term of $P(\alpha/2, y)$ where it is below the smallest normal double. Above $\alpha = 10^6$, $\gamma = 0$ is solved on the distribution function measured from the peak, as for any other tilt. At $\alpha = 1$ exactly the truncated-normal inverse is taken from `qnorm`, refined by one Newton step on the scaled complementary error function where the truncation point 0 lies three or more standard deviations above the mean $\gamma/(2\beta)$ of the untruncated normal, $-\gamma/\sqrt{2\beta} \ge 3$, and by Newton's method on that function alone where it lies more than 30 standard deviations above it. Near the truncation point it is solved by Halley's method on a short series, so that the quantile of a small lower-tail probability keeps its relative accuracy. ### A.5 Random generation at extreme parameters {-} The branches of the decision tree in §7 that start at extreme values come from the limits of double precision, and ordinary parameters never reach them. `?rmhn` states which sampler such a call uses and which calls stop with an error; this section gives the reasons. **Rounding in R's `rgamma`.** R's `rgamma`, from which the square-root-Gamma closed form and Sun et al.'s Algorithm 3 draw, rounds $\sqrt{a - 1/2}$ and the scale of a Gamma variate of shape $a$, and so shifts every variate alike by up to about $3 \times 2^{-53}\sqrt{a}$ of its standard deviation: $2.5 \times 10^{-7}$ at $a = 2^{59}$. Algorithm 3's own setup, which rounds the shape, the power and the rate of its Gamma proposal, adds a shift of the same kind; its draws were measured to shift by up to about $4.5 \times 2^{-53}\sqrt{\alpha}$ standard deviations in all ($3.8 \times 10^{-7}$ at $\alpha = 2^{59}$). The package keeps such common shifts to a few times $10^{-7}$ of a standard deviation. So the closed form, whose Gamma variate has shape $\alpha/2$, is used only up to $\alpha = 2^{60}$, and Algorithm 3 only up to $\alpha = 2^{59}$ (about $5.8 \times 10^{17}$). For the same reason, above $\alpha = 10^{18}$ Algorithm 1 is used with its Normal proposal whatever the constants of Theorem 1: its square-root-of-Gamma proposal would be built from $\sqrt{\alpha/2 - 1/2}$ rounded to double precision, a shift the acceptance test cannot see. Under a small tilt the Normal proposal accepts about $1/\sqrt{2}$ of its proposals there (0.706 at $\alpha = 2 \times 10^{18}$), below the 0.8 of Theorem 2e, and the law drawn stays exact. **Draws placed from the mode.** As $\alpha$ or the tilt grows, the law narrows around its mode relative to the mode itself, and a rounding error that moved every draw alike, such as that of a computed mode, would become a visible part of its standard deviation. Beyond the shapes above, the draws are therefore placed from the mode computed to about twice double precision, so that each is the exact draw rounded once, at any batch size: by the Normal proposal of Algorithm 1, centred at the mode, for a zero or negligible tilt of either sign ($|\gamma|/\sqrt{\beta}$ below `sqrt(.Machine$double.eps)`) above $\alpha = 2^{60}$, or by RTDR under `method = "rtdr"`; and by RTDR for any other $\gamma < 0$ above $\alpha = 2^{59}$. Sun et al. state Algorithm 1 for $\gamma > 0$ and Theorem 1 for $\gamma \geq 0$, but centred at the mode the Normal envelope holds for either sign. A positive tilt goes to Algorithm 1 above $\alpha = 2^{59}$, also at any batch size: no timing reaches that far, and RTDR cannot represent the law from a shape of about $10^{31}$ on (see below), while Algorithm 1 draws it. **Strong negative tilts.** Past a tilt $-\gamma/\sqrt{\beta} = 10^{150}$, where the law is its $\mathrm{Gamma}(\alpha, -\gamma)$ limit to double precision, Algorithm 3 proposes that limit itself and accepts practically every proposal, at a setup cost no envelope construction can undercut, so it is used there at any batch size. Its matching point for $\alpha \leq 1.1$ is then simply the mean $\alpha/|\gamma|$ of that limit, without a Newton step. **Where RTDR stops.** RTDR stops with an error for the whole call once the standard deviation of the law is no larger than the spacing of the doubles at its mode $m$, so that its draws could take only a few distinct values. Relative to $m$, the standard deviation is about $1/\sqrt{2\alpha}$ at $\gamma = 0$, and $1/\sqrt{\alpha}$ under a strong negative tilt, where the law is close to $\mathrm{Gamma}(\alpha, -\gamma)$; the doubles at $m$ are between $2^{-53} m$ and $2^{-52} m$ apart. So on the shape axis the stop comes between $2^{103}$ and $2^{105}$ ($4.1 \times 10^{31}$) at $\gamma = 0$, and between $2^{104}$ and $2^{106}$ ($8.1 \times 10^{31}$) under a strong negative tilt, depending on where $m$ falls between two powers of two. Under a positive tilt it comes from $\gamma/\sqrt{\beta}$ of about $9 \times 10^{15}$ ($2^{53}$ at $\beta = 1$). There the law is normal to far below double precision, with mean $m$ and variance $1/((\alpha - 1)/m^{2} + 2\beta)$, and can be drawn from directly. Algorithm 1 and the truncated normal for $\alpha \approx 1$ do not stop there: their draws are placed at the mode to within a few units in the last place, so that they take only a few distinct values near it. The default method uses Algorithm 1 for such laws with $\gamma > 0$ and $\alpha > 1$: up to $\alpha = 2^{59}$ their tilt exceeds $2\sqrt{\alpha}$, and above that shape Algorithm 1 serves at any tilt. **Where RTDR cannot place the law.** For $\alpha \le 1$ RTDR works with $\log(\sqrt{\beta} X)$, whose mode is the logarithm of about $\alpha\sqrt{\beta}/|\gamma|$ under a negative tilt. Below a shape of about $2.5 \times 10^{-324}\,|\gamma|/\sqrt{\beta}$, half the smallest positive double times the tilt, that number rounds to 0; this needs $\alpha$ below $4.4 \times 10^{-16}$, and almost all the mass of $X$ then lies below the smallest positive double. Algorithm 3 samples the law exactly there and is used instead, under `method = "rtdr"` too. RTDR cannot place the law either where $-\gamma/\sqrt{\beta}$ passes the largest double, or comes within a small factor of it at a shape of a few; `?rmhn` gives the sampler used, or the error raised, in each case. **Subnormal shapes.** `rmhn` draws the peak of a law with a subnormal shape (see `?rmhn`) through a proxy, the law of the shape $2^{-60}$ above $\sqrt{\beta}\,x = 1$, which has the peak's density there to within a relative $4 \times 10^{-18}$. ### A.6 Timings behind `method = "auto"` {-} §7 sets each sampler's setup against its cost per proposal. The implementation adds the following to those costs. - **Sun Algorithm 3.** For $\alpha > 1.1$ the initial matching point needs a bisection for the inflection point, and the single Newton refinement uses two `boost::math::digamma` evaluations and one `boost::math::trigamma`. The acceptance log-ratio reuses $T = (X/m)^{1/r}$ directly. Its three terms are each of the order of the Gamma rate $R$ and cancel to order one, so above $R = 10^{6}$, where their rounding would no longer be negligible, it is evaluated in the form $R\,w\,[(1+v)^{p} - 1 - p v]$, with $w = X/m$, $v = w - 1$ and $p = \beta m/(\beta m + |\gamma|)$. - **RTDR in region (a).** The search for each contact point starts from $m/2$ and $m + \max(1, m/2)$ — not from $m \mp \sqrt{-2/(\log f)''(m)}$, which Gao & Wang (2025) recommend after their Eq. (8) and use in the $\lambda \geq 1$ branch of the R code of their Appendix B — and walks in, in about $\log_2 m$ Newton steps, except where the construction is centred at the mode (once the magnitudes of the three terms of $\log f(m)$, which set the scale of its rounding error, sum past $0.01$ divided by the machine epsilon: from a shape of about $2.4 \times 10^{12}$ to $3.1 \times 10^{12}$ at moderate tilts, of about $2.3 \times 10^{13}$ near a tilt of $-\alpha/e$, and of about $4.5 \times 10^{13}$ at any tilt) or $\gamma/\sqrt{\beta}$ is below $-10^{10}\sqrt{\alpha}$; there it starts next to the mode, at the offset that Gao & Wang (2025, Section 2) recommend. Each proposal then takes three uniform variates and a few `log` evaluations. The crossovers come from two paired timing races of the two samplers, whose scripts and results are not shipped with the package, and from three runs of [inst/benchmarks/auto_dispatch.R](https://github.com/t-momozaki/mhn/blob/main/inst/benchmarks/auto_dispatch.R) (installed at `system.file("benchmarks", "auto_dispatch.R", package = "mhn")`), whose grid starts at $\alpha = 0.3$; all were timed on an Apple M1 Ultra with the code of this release, and the crossings move from machine to machine. The first race covers $\alpha$ from 0.01 to $10^4$ and $|\gamma|/\sqrt{\beta}$ from 0.5 to $10^4$; the second covers shapes from $10^3$ to $10^{14}$ under tilts of both signs, positive tilts from 1.25 to 3 times $\sqrt{\alpha\beta}$ at shapes from 1.5 to 1000, and single draws within one vectorised call, which pays a setup per draw without R's cost per call. For $\gamma < 0$ Algorithm 3 wins single draws for $\alpha < 1$ (by 11--23% at $\alpha = 0.01$, by at most 2% for $\alpha$ from 0.3 to 0.8), and RTDR takes over from $S = 5$, or $S = 25$ at $\alpha = 0.01$, where the spike at the origin makes Algorithm 3's setup relatively cheapest. For $\alpha > 1$ Algorithm 3 locates its matching point by bisection and Newton steps, and RTDR wins at every $S$ from $\alpha = 1.5$ to 100. From $\alpha = 10^3$, where RTDR's search walks in, its setup is the dearer: Algorithm 3 wins single draws, RTDR taking up to 1.4 times as long, and up to 5.2 times in a vectorised call, and RTDR takes over from $S = 5$ or 10 at $\alpha = 10^3$, 10 or 25 at $10^4$, 25 from $10^5$ to $10^7$, 25 or 50 at $10^8$ and 50 beyond. Where its search starts next to the mode, RTDR wins at every $S$, single draws included (by 3.6--5.0%, and by 17--34% in a vectorised call). For $\gamma > 0$, $\alpha > 1$ Algorithm 1 wins or ties single draws at every shape and tilt timed (RTDR takes 1.6 to 12 times as long in a vectorised call), and its per-proposal cost depends on $\gamma/\sqrt{\alpha\beta}$. Up to about 1 it proposes from the square root of a Gamma law, one `rgamma` call per proposal, and RTDR wins from $S = 5$ ($\alpha \le 5$), 10 ($\alpha = 10$), 25 ($\alpha = 100$), 25 or 50 ($\alpha$ from $10^3$ to $10^7$), 50 or 100 ($\alpha = 10^8$) or 50 to 200 ($\alpha = 10^{10}$ and $10^{12}$), by up to 69%. Above it, up to 2, Algorithm 1 proposes from a Normal law that accepts only 0.81--0.92 of its proposals, and RTDR still wins most large batches, from a larger $S$ and by up to 12% (32% at $\alpha = 10^7$, where Algorithm 1's acceptance test is dear); at the edge of that band, $\gamma/\sqrt{\alpha\beta} = 2$, it wins no $S$ at $\alpha = 1.5$, $10^{10}$ and $10^{12}$ (it is 3.0%, 2.0% and 2.3% slower at best). From about 2.5--3 on Algorithm 1 wins nearly every $S$: beyond 2, RTDR gains at most 2.3% at shapes up to 1000 and 0.4% at larger ones, but 12% at $\alpha = 10^7$. The rule cuts at $\gamma/\sqrt{\beta} = 2\sqrt{\alpha}$ and splits its thresholds at $\sqrt{\alpha}$, near Algorithm 1's own switch of proposal (at the tilt $\Delta_{\mathrm{crit}} \geq \sqrt{\alpha}$ where $K_{1} = K_{2}$, with $\Delta_{\mathrm{crit}}/\sqrt{\alpha} \to 1$ as $\alpha \to \infty$; Sun et al. 2023, Section 4.1.1): 5 and 25 below $\alpha = 10$, 25 and 50 up to $10^3$, 50 and 100 up to $10^8$, 100 and 1000 beyond, and 5 and 25 where RTDR's search starts next to the mode. ## Appendix B. Misprints in the source papers {-} ### B.1 Sun et al. (2023) {-} The published paper has six typesetting errors in results this vignette draws on that affect what is computed from them: four in the main text (two in Theorem 1, one in Lemma 7 and one in step 2 of Algorithm 3) and two in its Supplementary Material (the statement of Lemma 10, and the curvature $l''$ in Section 2.17). The paper has other slips in the results cited, which do not. In each of the six cases the paper itself carries the correct form elsewhere — in a proof, in a neighbouring display, or in the expression being differentiated — and the package uses the correct form. - **Square-root-of-Gamma density and selection rule in Theorem 1.** Theorem 1 gives the density of the square root of a Gamma variable as $f_{\sqrt{\mathrm{Gam}}}(x \mid \alpha, \delta) = x^{\alpha-1}\delta^{\alpha/2}e^{-\delta x^{2}}/\Gamma(\alpha/2)$, which integrates to $1/2$; the density is twice that, and the constant $K_{2}$ of its first part is the one for the density with the factor 2. The indicator that picks the proposal is printed as $I = \mathbb{I}(K_{1}(\delta, \alpha, \beta, \gamma) \le K_{2}(\mu, \alpha, \beta, \gamma))$, with the first arguments of $K_{1}$ and $K_{2}$ exchanged: $K_{1}$ is a function of $\mu$ and $K_{2}$ of $\delta$, as Algorithm 1 uses them. (The theorem also labels two of its parts "a".) The package draws the proposal as the square root of a Gamma variate, whose density carries the factor 2, and compares $K_{1}(\mu_{\mathrm{opt}}, \cdot)$ with $K_{2}(\delta_{\mathrm{opt}}, \cdot)$. - **Mixture weights** $p_{i}$ **in Lemma 7.** The main-text statement of Lemma 7 gives the weights for the $\sqrt{\mathrm{Gamma}}$ mixture representation as $p_{i} = \tfrac{1}{2\beta^{\alpha/2}} \cdot \tfrac{\Gamma((\alpha + i)/2)\, (\gamma/\sqrt{\beta})^{i}} {i!\, \Psi[\alpha/2,\, \gamma/\sqrt{\beta}]}$. By the definition of $\Psi$ these sum to $1/(2\beta^{\alpha/2})$ rather than one, contradicting the Poisson–Gamma hierarchical representation in Sun et al. (2023, Lemma 5a). The correct form, reproduced in the supplementary proof of Lemma 7 (Supplementary Section 2.13), is $$p_{i} \;=\; \frac{\Gamma((\alpha + i)/2)\, (\gamma/\sqrt{\beta})^{i}} {i!\, \Psi[\alpha/2,\, \gamma/\sqrt{\beta}]}$$ i.e.\ drop the leading $1/(2\beta^{\alpha/2})$. This affects Algorithm 2 only, which the package does not implement; the correction is recorded here for completeness. - **Series truncation discriminants** $C_{1}$, $C_{2}$ **in Lemma 10.** Lemma 10 (Supplementary Section 1) splits the series for $\Psi[\alpha/2, x]$ into its even and odd terms, $A(k)$ and $B(k)$. Its parts (a) and (b) give indices $C_{1}$ and $C_{2}$ beyond which the ratio of successive terms, $A(k+1)/A(k)$ or $B(k+1)/B(k)$, stays below a fixed contraction factor $q \in (0, 1)$, and part (d) sets the truncation point $K = \max\{K_{1}, K_{2}\}$ for a tolerance $\epsilon$ at the indices, no smaller than $C_{1}$ and $C_{2}$, past which every term is at most $\epsilon/4$; the package truncates at the upper bound on them that the geometric decay gives. Reducing the ratio bound to a quadratic in $k$ gives a discriminant whose last term is $\alpha x^{2}$ (and $(\alpha+1) x^{2}$ for $C_{2}$). The statement of Lemma 10 prints these as $\alpha x$ and $(\alpha+1) x$ — a single missing power of $x$. Its proof (Supplementary Sections 2.2.1 and 2.2.2) carries $z^{2}$ throughout, in agreement with the rederivation. The corrected $x^{2}$ form, with $q = 1/2$, is what the package uses in [src/mhn_psi_series.cpp](https://github.com/t-momozaki/mhn/blob/main/src/mhn_psi_series.cpp) (marked with an `// ERRATA` comment). With the printed form, $C_{1}$ and $C_{2}$ can fall below the index from which the ratio is actually at most $q$, so the geometric tail bound of part (d) no longer follows. The package uses the corrected form, which is the one the proof supports. - **Gamma shape in step 2 of Algorithm 3.** Algorithm 3 proposes $X = m T^{r}$ with $r = (\beta m + |\gamma|)/(2\beta m + |\gamma|)$. The display that introduces it gives the Gamma shape of $T$ as $\alpha(\beta m + |\gamma|)/(2\beta m + |\gamma|) = \alpha r$, which makes $X$ follow the proposal kernel of Theorem 3; step 2 prints it as $(\beta m + |\gamma|\alpha)/(2\beta m + |\gamma|)$, with the parentheses lost. The two agree only at $\alpha = 1$, and with the printed shape the acceptance test of step 5 would no longer return MHN draws. The package uses $\alpha r$. The same display sets $Y = T^{r}$, whose density matches the kernel of Theorem 3 only at $m = 1$; step 3's $X = m T^{r}$, which the package draws, matches it for every $m$. - **Sign in the curvature $l''$ of Supplementary Section 2.17.** The Newton–Raphson refinement of the Algorithm 3 matching point uses the first two derivatives of the log acceptance probability $l(m)$. Differentiating the printed $l'$ gives a first term $-4\alpha\beta^{2}|\gamma|/(2\beta m + |\gamma|)^{3}$ times $\psi(\alpha(\beta m + |\gamma|)/(2\beta m + |\gamma|)) - \log(\beta m^{2} + m|\gamma|)$; the printed $l''$ has it with a plus sign. At the starting point $m_{\mathrm{init}}$ that changes the length of the Newton step: at $(\alpha, \beta, \gamma) = (3, 1, -2)$, with $m_{\mathrm{init}} = 0.772$, $l''$ is $-1.03$ and the printed form gives $-1.17$. Farther from it the printed curvature can take the wrong sign, which reverses the step: at $(2, 1, -0.1)$ and $m = 0.223$, $l''$ is $-20.9$ and the printed form gives $+1.9$. The package uses the minus sign. ### B.2 Gao & Wang (2025) {-} Their Eq. (15), the envelope of region (d), prints the left-tangent interval as $(y^{*}, p_l)$, leaving $(\dot t_l, y^{*}]$ uncovered, where $\dot t_l < y^{*}$ is the point at which $\log g$ has the slope of the tangent at $t_l$; the package follows their Theorem 4.3, $(\dot t_l, p_l]$, and like their Appendix B code starts the tangent at the last breakpoint kept.