Pith. sign in

REVIEW 3 major objections 4 minor 14 references

Complexity of exact sampling of the first passage of a stable subordinator

T0 review · 3 major / 4 minor · reviewed 2026-08-07 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗

Pith's one-line read This paper proves that exact sampling of the first passage of a stable subordinator can be done with expected complexity $O(1+|\ln(1-\alpha)|)$, using acceptance-rejection envelopes rather than numerical inversion or integration.

desk verdict The chi sampler is a real algorithmic advance, but the abstract's 'without numerical inversions' claim for the full first-passage algorithm quietly assumes an oracle for B^{-1}. read the letter →

arxiv 2506.03047 v1 pith:BOSGS7SE submitted 2025-06-03 stat.CO

Zhiyi Chi This is my paper · ORCID
classification stat.CO MSC 65C0560G5160G5268Q25
keywords exactsamplingfirstpassagestablesubordinatoracceptance-rejectioncomplexityLévyprocessbarriergammadistribution
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

The paper claims that the first passage of a stable subordinator across a non-increasing regular barrier can be sampled exactly with expected complexity $O(1+|\ln(1-\alpha)|)$, where $\alpha$ is the index of the subordinator. It reduces the sampling problem to drawing from a bivariate distribution parameterized by $\alpha$ and a scalar $z$ that is independent of the barrier, then gives three acceptance-rejection algorithms covering different regions of $(\alpha,z)$. The hardest region, $\alpha$ close to 1 and $z$ close to 0, is handled by an algorithm whose expected running time is $O(|\ln(1-\alpha)|)$, built from piecewise envelopes. If the claim is right, exact simulation of the first passage no longer needs numerical inversion or integration, and the earlier $\delta^{-3}$ blow-up near $\alpha=1$ is replaced by a logarithmic cost.

What carries the argument

The central object is the normalized density $\chi_{\alpha,z}$ and the ratio function $H_\alpha(\theta)$ defined through sinc functions, which is strictly increasing, satisfies the ratio-symmetry property that $\alpha^{-1}\ln H_\alpha(\theta)$ is symmetric about $\alpha=1/2$, and has monotonicity bounds developed in Proposition 2.2. The paper builds a piecewise-constant approximation $J_\alpha(\theta)$ of $H_\alpha(\theta)$ on $[0,\pi)$ with $O(\ln(1/\delta))$ pieces, each piece within factor $1+\Delta$ of $H_\alpha$, so that the envelope $Q_{\alpha,z}$ becomes easy to integrate and sample. The elementary building blocks are explicit envelopes $\varphi^*_{a,b,c}$ and $g^*_{a,b,c}$ for functions $x^{c-1}e^{\pm x}$ on intervals $(a,b]$, and the logarithmic map $v=\ln(1+y)$ converts the $y$-part of $\chi_{\alpha,z}$ into mixtures of these elementary functions, allowing exact sampling through acceptance-rejection.

What would settle it

A concrete test: implement the combined algorithm and measure the expected number of iterations for $\alpha=1-10^{-k}$ with $k=1,\ldots,8$ and $z$ as small as $10^{-60}$; the claim predicts growth at most linear in $k$, so observing growth like $2^k$ or failure of the Gamma sampler with shape $\delta=10^{-k}$ to terminate in bounded time would falsify the central complexity claim.

Watch

Extended reading notes

Core claim

The central claim is that exact sampling of the first passage reduces to sampling the normalized bivariate density $\chi_{\alpha,z}(y,\theta) = [1-(y+1)^{-\delta/\alpha}]^{-\alpha} H_\alpha(\theta) \exp(-z H_\alpha(\theta)(y+1))$, where $\delta=1-\alpha$ and $H_\alpha(\theta)=\mathrm{sinc}(\delta\theta)(\mathrm{sinc}(\alpha\theta)/\mathrm{sinc}(\theta))^{\alpha/\delta}$, with $z$ a barrier-independent random variable represented as $z=\xi/H_\alpha(\Theta)$ for $\xi\sim\mathrm{Exp}(1)$ and $\Theta\sim\mathrm{Unif}(0,\pi)$. The paper establishes that the normalized $\chi_{\alpha,z}$ can be sampled with expected complexity $O(1+|\ln\delta|)$ uniformly in $z>0$, using acceptance-rejection with explicit envelopes and no numerical inversion or integration. Three algorithms are combined: one with complexity $O(z^{-1}(1+\delta/z)^2)$, one with complexity $O(\delta^{-2}/C_z)$ with $C_z>0$ increasing in $z$, and a third working near $(1,0)$ that samples the normalized function $Q_{\alpha,z}(\theta)=\psi_\alpha(zH_\alpha(\theta))e^{-zH_\alpha(\theta)}$ via a piecewise-constant proxy $J_\alpha(\theta)$ for $H_\alpha(\theta)$ accurate within factor $1+\Delta$ on $O(\ln(1/\delta))$ intervals.

Load-bearing premise

The $O(1+|\ln\delta|)$ bound assumes that the building-block samplers—Gamma distributions with shape $\delta\to 0$ and the mixture pieces in Section 4—have uniformly bounded expected running time, which the paper cites from [6,11] without proof; it also assumes $B^{-1}(s)$ in Algorithm 2.1 is available without numerical inversion, a point not covered by the abstract's claim of no numerical inversions.

Editorial extensions

If this is right

  • The first passage of a stable subordinator can be sampled exactly with $O(1+|\ln(1-\alpha)|)$ expected operations measured in uniform random bits and closed-form function evaluations, for non-increasing regular barriers.
  • Sampling $S_{\tau-}$ conditioning on $\tau$ and passage by jump is reduced to a single barrier-independent bivariate sampling problem, so changing the barrier only requires evaluating $B^{-1}(s)$.
  • The representation $z=\xi/H_\alpha(\Theta)$ shows that difficulty concentrates near $\alpha\to 1$ and $z\to 0$, and the complexity degrades only logarithmically in that region.
  • The algorithm produces samples from the same distribution as the earlier fast algorithm, and in experiments it is dramatically faster for ultra-small $z$ such as $z\le 10^{-50}$.
  • Because no numerical inversion or integration is used, the complexity does not scale with the number of precision bits $N$ required of the output, under the paper's idealization of closed-form functions.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • The same piecewise-constant envelope construction of $J_\alpha$ may transfer to related first-passage or hitting problems for tempered stable subordinators, where a similar $H$-type ratio function would need an analogous approximation.
  • The uniform-complexity assumption for small-shape Gamma samplers could be tested in isolation by benchmarking expected running time as the shape parameter $\delta$ tends to $0$; if the bound fails, the overall $O(1+|\ln\delta|)$ claim would degrade.
  • A concrete stress test would run the sampler for $\alpha=1-10^{-k}$ with $k$ up to about 12 and for $z$ as small as $10^{-60}$, and check that the expected number of iterations grows linearly in $k$ rather than polynomially or exponentially in $1/\delta$.
  • The paper leaves implicit that the same reduction to $\chi_{\alpha,z}$ would apply to barriers that are only piecewise smooth, provided the regularity conditions on $b$ in Section 2 are met locally.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

3 major / 4 minor

Summary. The paper studies exact sampling of the first passage of a stable subordinator of index α across a non-increasing regular barrier. It reduces the problem to sampling a normalized two-dimensional function χ_{α,z}(y,θ), where z is a barrier-independent scalar, and develops three acceptance-rejection algorithms. The first two (§3) are simple and have complexity O(1+z^{-1}(1+δ/z)^2) and O(δ^{-2}/C_z). The third (§§4–6) uses a logarithmic transform y=ℓ(v), piecewise-constant approximations of H_α, and envelope constructions for log-concave and Gamma-restricted densities to obtain expected complexity O(|ln δ|) for (α,z) close to (1,0). Combining the algorithms gives the headline bound O(1+|ln(1-α)|) for sampling from the normalized χ_{α,z}, and hence, via Algorithm 2.1, for exact sampling of the first passage. The paper includes detailed proofs in the Appendix, numerical comparisons, and R code on GitHub.

Significance. If the stated gaps are closed, this is a substantial improvement over the earlier algorithm of González Cázares, Lin, and Mijatović, which had complexity O(δ^{-3}+|ln α|+ln N) and relied on numerical inversion and integration. The reduction to the barrier-independent quantity z and the normalized χ_{α,z} is clean, and the envelope constructions in Sections 4–6 are original and technically involved. The paper also provides reproducible code and empirical comparisons among the three algorithms. The main technical result on sampling χ_{α,z} appears sound and is a genuine contribution. However, the headline claim about first-passage sampling without numerical inversions is not fully established as stated, because Algorithm 2.1 requires evaluating B^{-1}, and the paper does not provide B^{-1} as a closed-form oracle for general barriers.

major comments (3)
  1. [Abstract, §1, and Algorithm 2.1 (line 2)] The paper claims exact sampling of the first passage with complexity O(1+|ln δ|) 'without numerical inversions or integrations'. However, Algorithm 2.1 requires t ← B^{-1}(s), where B(t)=t^{-1/α}b(t) is defined in (2.2). For a general non-increasing regular barrier b, B^{-1} is not a closed-form function; it would require numerical root-finding, which is outside the exact-sampling model stated in Section 1. The experiments in Section 7.3 use only two barriers with explicit inverse, B^{-1}(s)=(10/s)^α and B^{-1}(s)=[100/(s+1)]^α, so they do not test the general claim. The 'without numerical inversions' assertion should be restricted to the χ_{α,z} sampler, or the paper should state that B^{-1} is provided as an exact oracle and count that oracle call in the complexity.
  2. [§4.1, equation (4.9)] The proof of (4.9) relies on the assertion that r_n(4) ≤ 2(n+1)/(n+2) for n≥7, justified only by 'It is not hard to show, either analytically or numerically'. This inequality is needed to ensure ∫ f_n ≤ 2M(a,A,c), which enters the universal ratio 4 in (4.13) and hence the correctness of the envelope φ*_{a,b,c}. A load-bearing step is left unproved; please supply a proof or a precise reference.
  3. [§3.1–3.2 and §7.1] The uniform complexity bounds in Algorithms 3.1 and 3.2 assume that sampling from a Gamma density has uniformly bounded expected time for all shape parameters, including δ=1-α arbitrarily close to 0. Sections 3.1 and 3.2 cite [6] for this uniformity, and Section 7.1 notes that numerical issues for small δ are handled by sampling the logarithm of the Gamma variate. The O(1+|ln δ|) claim depends on this uniformity, so the paper should state the precise theorem or algorithm in [6] that guarantees it, and confirm that it applies within the paper's exact-sampling model rather than in an approximate numerical sense.
minor comments (4)
  1. [Proposition 2.1] The statement 'For θ∈(0,π) and θ∈(0,1)' contains a typo: the second interval should be for α, not θ.
  2. [Algorithms 4.1 and 4.2] The instructions 'Sample Y∼E_s' and 'Sample Z∼E_s' use E_s, which is defined in (4.1) as a density, not as a distribution name. Please specify the distribution associated with the density E_s.
  3. [Algorithm 7.1] The pseudocode returns (z,t,y,θ), but unlike Algorithm 2.1 it does not return the computed undershoot x or the jump variable. The text explains this, but a brief comment in the pseudocode would make the output format clearer.
  4. [§7.1] The notation 'ln #γ_{δ,τ}' is introduced without an explicit definition; since the pushforward notation was defined in Section 4, defining it here as that pushforward under the logarithmic map would improve readability.

Circularity Check

0 steps flagged · score 2.0 of 10

No material circularity; the central O(1+|ln(1−alpha)|) sampling bound is derived from explicit envelope inequalities, with self-citations used only as external prior results.

full rationale

The paper's central claim is a constructive complexity bound for sampling normalized chi_{alpha,z} by acceptance-rejection. Every step in Sections 3, 5 and 6 gives explicit envelopes (A_{alpha,z}, B_{alpha,z}, G_{alpha,z}, P_{alpha,z}, Q_{alpha,z}) and derives acceptance-probability lower bounds from the paper's own Proposition 2.1/2.2 and sinc-series identities; no fitted parameter is renamed as a prediction. The reduction from first-passage sampling to chi_{alpha,z} is algebraic (Section 2, equations (2.2)-(2.5)) and uses [3] and [8] only for the known law tau ~ B^{-1}(zeta) and the conditional undershoot density; these are prior published results, not results of the present paper being invoked as their own support. The author's self-citation [3] is thus minor and not load-bearing in a circular sense. The O(|ln delta|) bound does depend on uniform boundedness of Gamma and elementary building-block samplers cited from [6,11]; this is an external-support assumption and, if false, would weaken the theorem but does not make the derivation circular. One non-circular limitation should be noted: Algorithm 2.1 still calls B^{-1}(s), so for a general barrier the 'without numerical inversions' headline is not fully established; this is an incompleteness/correctness concern, not a circularity.

Assumptions & free parameters 0 free parameters · 5 assumptions · 0 invented entities

No free parameters fitted to data. The constants Δ, α_0, θ_0 in (6.1) are fixed arbitrary choices within given ranges; the O(|ln δ|) bound holds uniformly with constants depending on them. The axioms are standard published results on stable subordinators and sampling subroutines.

assumptions (5)
  • domain assumption First passage time representation τ ~ B^{-1}(ς) with ς having stable density p (Eq. 2.2).
    Taken from Chi (2016) [3]; not proved in this paper, used as the starting reduction.
  • domain assumption Conditional distribution of undershoot and jump at passage (Eq. 2.3), including creeping probability, from [3,8].
    Provides the target density χ_{α,z}.
  • standard math Stable density integral representation (2.1) from Sato [13] and Zolotarev [14].
    Used to derive H_α and χ.
  • domain assumption Uniformly bounded exact samplers for Gamma and normal distributions exist (Devroye [6], Liu et al. [11]).
    Used inside AR loops; the complexity bound inherits this.
  • ad hoc to paper Closed-form functions (power, exp, log, trig, Gamma) are evaluated exactly in O(1) time.
    This is the paper's explicit conceptualization of exact sampling, stated in Section 1.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Complexity of exact sampling of the first passage of a stable subordinator." pith.science (2026). https://pith.science/paper/BOSGS7SE

@misc{pith2026250603047,
  author       = {Pith},
  title        = {Pith review of: Complexity of exact sampling of the first passage of a stable subordinator},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/BOSGS7SE}},
  note         = {Machine review of arXiv:2506.03047}
}
abstract

We consider the exact sampling of the first passage of a stable subordinator across a non-increasing regular barrier. First, the sampling is reduced to one from a bivariate distribution parameterized by the index $\alpha$ of the subordinator and a scalar $z$ independent of the barrier. Then three algorithms are devised for different regions of $(\alpha, z)$, using the acceptance-rejection method without numerical inversion or integration. When combined, the algorithms allow the exact sampling of the first passage to be done with complexity $O(1+|\ln(1-\alpha)|)$.

Figures

Figures reproduced from arXiv: 2506.03047 by the authors.

Figure 1
Figure 1. Comparison of Algorithms 3.1 and 3.2 when [PITH_FULL_IMAGE:figures/full_fig_p025_1.png] view at source ↗
Figure 2
Figure 2. Comparison of Algorithms 3.1 and 3.2 under the same setting as in Figure 1 [PITH_FULL_IMAGE:figures/full_fig_p025_2.png] view at source ↗
Figure 3
Figure 3. Comparison of running times of Algorithms 3.1 and 3.2 to sample 2 [PITH_FULL_IMAGE:figures/full_fig_p027_3.png] view at source ↗
Figures from the paper (3 more)
Figure 4
Figure 4. Figure 4: Comparison of Algorithms 3.2 and 5.1 when [PITH_FULL_IMAGE:figures/full_fig_p027_4.png]
Figure 5
Figure 5. Figure 5: Comparison of Algorithms 3.2 and 5.1 under the same setting as in Figure 4 [PITH_FULL_IMAGE:figures/full_fig_p029_5.png]
Figure 6
Figure 6. Figure 6: Mean and SD of running time (in s) of Algorithm 7.1 over 100 repetitions, with [PITH_FULL_IMAGE:figures/full_fig_p029_6.png]

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

14 extracted references · 14 canonical work pages

  1. [8]

    Fast exact simulation of the first passage of a tempered stable subordinator across a non-increasing function

    Gonz \'a lez C \'a zares, J. , Lin, F. , and Mijatovi \'c , A. (2023). Fast exact simulation of the first passage of a tempered stable subordinator across a non-increasing function. arXiv\/ , arXiv:2303.11964

    work page Pith review arXiv 2023

  2. [6]

    Devroye, L. (1986). Nonuniform Random Variate Generation . Springer-Verlag, New York

    work page 1986

  3. [1]

    Asmussen, S. r. and Albrecher, H. (2010). Ruin probabilities\/ , Second ed. Advanced Series on Statistical Science & Applied Probability, Vol. 14. World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ

    work page 2010

  4. [2]

    Chambers, J. M. , Mallows, C. L. , and Stuck, B. W. (1976). A method for simulating stable random variables. J. Amer. Statist. Assoc.\/ 71, 354, 340--344

    work page 1976

  5. [3]

    Chi, Z. (2016). On exact sampling of the first passage event of a L \'evy process with infinite L \'evy measure and bounded variation. Stochastic Processes Appl.\/ 126, 4, 1124--1144

    work page 2016

  6. [4]

    Chi, Z. (2025). Sampling first passage of stable subordinator. R package at https://github.com/zhc05001/Sampling-first-passage-of-stable-subordinator

    work page 2025

  7. [5]

    , Lim, J

    Dassios, A. , Lim, J. W. , and Qu, Y. (2020). Exact simulation of a truncated L \' e vy subordinator. ACM Trans. Model. Comput. Simul.\/ 30, 3, Art. 17, 26

    work page 2020

  8. [7]

    Devroye, L. (2012). A note on generating random variables with log-concave densities. Statist. Probab. Lett.\/ 82, 5, 1035--1039

    work page 2012

Show all 14 references
  1. [9]

    , Leydold, J

    H \"o rmann, W. , Leydold, J. , and Derflinger, G. (2004). Automatic Nonuniform Random Variate Generation . Statistics and Computing. Springer-Verlag, Berlin

    2004

  2. [10]

    Kyprianou, A. E. (2014). Fluctuations of L \'evy Processes with Applications\/ , Second ed. Universitext. Springer, Heidelberg

    2014

  3. [11]

    , Martin, R

    Liu, C. , Martin, R. , and Syring, N. (2017). Efficient simulation from a gamma distribution with small shape parameter. Computational Statistics\/ 32 , 1767--1775

    2017

  4. [12]

    Olver, F. W. J. , Lozier, D. W. , Boisvert, R. F. , and Clark, C. W. , Eds. (2010). NIST Handbook of Mathematical Functions . U.S. Department of Commerce National Institute of Standards and Technology, Washington, DC

    2010

  5. [13]

    Sato, K.-I. (1999). L\'evy Processes and Infinitely Divisible Distributions . Cambridge Studies in Advanced Mathematics, Vol. 68. Cambridge University Press, Cambridge

    1999

  6. [14]

    Zolotarev, V. M. (1966). On the representation of stable laws by integrals. Selected Translations in Mathematical Statistics and Probability\/ 6 , 84--88. First published in Russian in 1964

    1966

Pith tools

Reviewed August 7, 2026 · model on record in the stance chip above.