pith. sign in

arxiv: 2601.19633 · v3 · submitted 2026-01-27 · 🧮 math.PR · cs.NA· math.NA

Computing the density of the Kesten-Stigum limit in supercritical Galton-Watson processes

Alice Cortinovis , Sophie Hautphenne , Stefano Massei This is my paper

Pith reviewed 2026-05-16 10:55 UTC · model grok-4.3

classification 🧮 math.PR cs.NAmath.NA
keywords Galton-Watson processKesten-Stigum limitLaplace-Stieltjes transformLaguerre polynomialsmoment matchingdensity approximationbranching processessupercritical processes
0
0 comments X

The pith

A functional equation for the Laplace-Stieltjes transform combined with Laguerre polynomial moment matching yields accurate densities for the Kesten-Stigum limit in supercritical Galton-Watson processes.

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

This paper develops a numerical technique to compute the probability density of the limit random variable that arises in supercritical Galton-Watson branching processes. The method begins from the functional equation satisfied by the Laplace-Stieltjes transform of the limit law and approximates the density as a finite linear combination of Laguerre polynomials damped by an exponential factor. Coefficients are selected so the first several moments of the approximation match those implied by the functional equation. The procedure is tested on cases where the offspring generating function is a polynomial of bounded degree. The resulting densities quantify the random scaling factor that multiplies the deterministic exponential growth of the population.

Core claim

The density of the Kesten-Stigum limit can be approximated by solving the functional equation for its Laplace-Stieltjes transform through a moment-matching procedure that expresses the density as a linear combination of Laguerre polynomials multiplied by an exponential damping factor. This yields accurate numerical estimates when the offspring distribution has a polynomial generating function.

What carries the argument

The functional equation characterizing the Laplace-Stieltjes transform of the Kesten-Stigum limit law, solved via moment matching within the span of exponentially damped Laguerre polynomials.

If this is right

  • The density becomes computable for any supercritical Galton-Watson process whose offspring generating function is a polynomial of bounded degree.
  • Higher numbers of Laguerre terms produce successively better approximations by incorporating additional moments from the functional equation.
  • The method supplies a deterministic alternative to direct simulation for studying the distribution of the random growth amplitude.
  • Validation on several polynomial examples confirms that the approximation reproduces known properties of the limit law in those cases.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • The same moment-matching strategy might extend to offspring distributions with infinite support provided the functional equation can generate the required moments.
  • Analogous transform-based approximations could apply to other limit laws in branching processes that obey functional equations for their Laplace transforms.
  • Numerical tests on non-polynomial cases would reveal whether the exponential damping parameter requires adjustment for stability.
  • The approach could be compared against other basis expansions to determine the relative efficiency of Laguerre polynomials for this problem.

Load-bearing premise

The Laguerre polynomial expansion with moment matching produces stable and accurate density estimates for the limit distribution when the offspring generating function is a polynomial of bounded degree.

What would settle it

Generate a large Monte Carlo histogram of the normalized population size after many generations for a Poisson offspring distribution with mean greater than one and check whether the Laguerre approximation converges to the histogram shape as the number of terms increases.

Figures

Figures reproduced from arXiv: 2601.19633 by Alice Cortinovis, Sophie Hautphenne, Stefano Massei.

Figure 1
Figure 1. Figure 1: Relative residual Res(φ (k) ) along the iterations of the fixed-point method and Newton’s method when solving (4) for P(z) = 1 10 (1 + z + z 2 + · · · + z 9 ). see Remark 4. As a result, Newton’s method is faster than the fixed-point iteration by a factor between 10 and 30 when m ∈ {1.1, 1.25}. For larger values of m, the execution times of the two methods become comparable. The execution time of the forwa… view at source ↗
Figure 2
Figure 2. Figure 2: Condition numbers as a function of the basis dimension [PITH_FULL_IMAGE:figures/full_fig_p018_2.png] view at source ↗
Figure 3
Figure 3. Figure 3: Approximations of the density of W for GW processes with offspring p.g.f. P1(z) (left) and P2(z) (right). In both cases, the extinction probability is q = 0, and hence there is no point mass at 0. 19 [PITH_FULL_IMAGE:figures/full_fig_p019_3.png] view at source ↗
Figure 4
Figure 4. Figure 4: Approximations of the density of W for GW processes with offspring p.g.f. P3(z) (left) and P4(z) (right). For P3(z), a smaller number of generations (T = 8) was used in the simulation. In both cases, the extinction probability is strictly positive; the red curve represents our approximation of the absolutely continuous part of the density. 4.3 Comparison with generalized Gamma distribution fitting Another … view at source ↗
Figure 5
Figure 5. Figure 5: Comparison of the Laguerre-based approximation ( [PITH_FULL_IMAGE:figures/full_fig_p021_5.png] view at source ↗
Figure 6
Figure 6. Figure 6: Plot corresponding to the density of W for the whooping crane. For this example, we have m ≈ 1.04 and the distribution has a heavier right tail compared to the previous examples. For this reason, to estimate the parameter β we have used the 40%-to-70% range of the histogram, resulting in β ≈ 0.246. The simulation was run until time T = 100. 4.5 Handling non-polynomial offspring probability generating funct… view at source ↗
Figure 7
Figure 7. Figure 7: Plot corresponding to the density of W for the Chatham Islands black robin population. For this example, we have m ≈ 1.68 and the estimated parameter β is 0.543. The simulation was run until T = 15. approximation. A possible explanation for this behavior is that the coefficients of the offspring generating function decay more quickly when c is small, making the truncated generating function “closer” to a p… view at source ↗
Figure 8
Figure 8. Figure 8: Approximate density of the time at which the population size first reaches [PITH_FULL_IMAGE:figures/full_fig_p024_8.png] view at source ↗
Figure 9
Figure 9. Figure 9: Approximations of the coefficients of φ obtained by Newton’s method and the forward substi￾tution method for the two linear fractional examples considered in Section 4.5. Note that the forward substitution method achieves very precise approximations in both cases. [11] A. Edelman and G. Strang. Pascal matrices. The American Mathematical Monthly, 111(3):189–197, 2004. [12] J. Fernley and E. Jacob. A univers… view at source ↗
read the original abstract

This paper proposes a novel numerical method for computing the density of the limit random variable associated with a supercritical Galton-Watson process. This random variable captures the effect of early demographic fluctuations and determines the random amplitude of long-term exponential population growth. While the existence of a non-trivial limit is ensured by the Kesten-Stigum theorem, computing its density in a stable and efficient manner for arbitrary offspring laws remains a significant challenge. The proposed approach leverages a functional equation that characterizes the Laplace-Stieltjes transform of the limit distribution and combines it with a moment-matching method to obtain accurate approximations within a class of linear combinations of Laguerre polynomials with exponential damping. The effectiveness of the approach is validated on several examples in which the offspring generating function is a polynomial of bounded degree.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

2 major / 1 minor

Summary. The paper proposes a numerical method to compute the density of the Kesten-Stigum limit random variable W for supercritical Galton-Watson processes. It solves the functional equation satisfied by the Laplace-Stieltjes transform of W and recovers an approximation to the density by moment-matching within the span of Laguerre polynomials multiplied by an exponential damping factor. The approach is illustrated and validated only on examples where the offspring probability generating function is a polynomial of bounded degree.

Significance. If the method can be shown to be stable and accurate for general offspring distributions (including non-polynomial cases such as Poisson), it would supply a practical computational tool for the distribution of the limiting multiplier W, which governs the random amplitude of exponential growth in branching processes and appears in many applied models. The construction is non-circular and rests directly on the classical Kesten-Stigum functional equation.

major comments (2)
  1. Numerical experiments section: validation is confined to polynomial offspring laws of bounded degree. For general laws the higher moments must be obtained recursively from the functional equation, yet no results, error bounds, or stability analysis are supplied for non-polynomial cases (e.g., Poisson or geometric offspring). This directly limits support for the claim of applicability to arbitrary offspring distributions.
  2. Laguerre approximation paragraph: the truncation degree of the Laguerre expansion and the damping parameter are free parameters whose selection is not accompanied by any convergence-rate analysis or a priori error estimate. When the density of W is less regular than in the polynomial-offspring case, the rate at which the series converges is left unexamined.
minor comments (1)
  1. Abstract and introduction: the phrase 'arbitrary offspring laws' should be qualified to reflect that all reported numerical evidence is restricted to polynomial generating functions.

Simulated Author's Rebuttal

2 responses · 0 unresolved

We thank the referee for the careful reading and constructive suggestions. We address the two major comments point by point below and indicate the revisions we will incorporate.

read point-by-point responses

  1. Referee: Numerical experiments section: validation is confined to polynomial offspring laws of bounded degree. For general laws the higher moments must be obtained recursively from the functional equation, yet no results, error bounds, or stability analysis are supplied for non-polynomial cases (e.g., Poisson or geometric offspring). This directly limits support for the claim of applicability to arbitrary offspring distributions.

    Authors: We agree that the present numerical section is restricted to polynomial offspring distributions. The method itself is formulated for arbitrary offspring laws: the functional equation holds generally, and all moments of W can be obtained recursively from the offspring distribution whenever the Kesten-Stigum condition is satisfied. In the revised manuscript we will add a new subsection containing numerical results for the Poisson offspring distribution (mean 2), including a comparison of the recovered density against the known mean and variance of W, together with a brief stability check on the moment recursion. This will directly address the concern about general applicability. revision: yes

  2. Referee: Laguerre approximation paragraph: the truncation degree of the Laguerre expansion and the damping parameter are free parameters whose selection is not accompanied by any convergence-rate analysis or a priori error estimate. When the density of W is less regular than in the polynomial-offspring case, the rate at which the series converges is left unexamined.

    Authors: We acknowledge that a rigorous a priori error analysis is currently missing. The truncation degree is chosen in practice by monitoring the residual of the functional equation, and the damping parameter is fixed by matching the known exponential decay rate of the density of W. For polynomial offspring distributions the density is C^∞, producing spectral convergence; for general laws the rate is limited by the smoothness of the density. In the revision we will add a short discussion of the expected convergence behavior in terms of Sobolev regularity of the target density and will report empirical convergence rates (residual versus degree) for both the existing polynomial examples and the new Poisson case. revision: yes

Circularity Check

0 steps flagged

No significant circularity; derivation rests on established Kesten-Stigum equation plus independent Laguerre moment-matching

full rationale

The paper invokes the standard Kesten-Stigum functional equation for the Laplace-Stieltjes transform of the limit W (a classical result, not derived here) and computes moments from it before applying a standard Laguerre-polynomial moment-matching approximation with exponential damping. Neither step reduces the claimed density to a fitted parameter or to a self-referential definition; the approximation class is chosen independently of the target density. Validation is restricted to polynomial offspring pgfs, but this is a limitation on demonstrated accuracy, not a circularity in the derivation chain. No self-citation is load-bearing for the central claim, and the method remains self-contained against external benchmarks.

Axiom & Free-Parameter Ledger

1 free parameters · 1 axioms · 0 invented entities

The approach rests on the theoretical functional equation from prior literature and introduces a numerical approximation scheme with one main free parameter for the polynomial degree.

free parameters (1)
  • truncation degree for Laguerre expansion
    Parameter controlling approximation accuracy, chosen based on validation examples.
axioms (1)
  • domain assumption The limit random variable satisfies the functional equation from the Kesten-Stigum theorem
    Invoked as the basis for the Laplace-Stieltjes transform characterization.

pith-pipeline@v0.9.0 · 5438 in / 1104 out tokens · 22494 ms · 2026-05-16T10:55:36.183135+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

20 extracted references · 20 canonical work pages

  1. [1]

    K. B. Athreya and P. E. Ney. Branching processes, volume 196. Springer Science & Business Media, 1972

    work page 1972

  2. [2]

    Baker, P

    J. Baker, P. Chigansky, K. Hamza, and F. C. Klebaner. Persistence of small noise and random initial conditions. Advances in Applied Probability , 50(A):67–81, 2018

    work page 2018

  3. [3]

    A. D. Barbour, P. Chigansky, and F. C. Klebaner. On the emergence of random initial conditions in fluid limits. Journal of Applied Probability , 53(4):1193–1205, 2016

    work page 2016

  4. [4]

    Bauman, P

    N. Bauman, P. Chigansky, and F. Klebaner. An approximation of populations on a habitat with large carrying capacity. arxiv. arXiv preprint arxiv:2303.03735 , 2023

    work page arXiv 2023

  5. [5]

    Braunsteins, S

    P. Braunsteins, S. Hautphenne, and J. Kerlidis. Existence and non-existence of consistent estimators in supercritical controlled branching processes. arXiv preprint arXiv:2504.03389 , 2025

    work page arXiv 2025

  6. [6]

    Braunsteins, S

    P. Braunsteins, S. Hautphenne, and C. Minuesa. Consistent least squares estimation in population- size-dependent branching processes. Journal of the American Statistical Association , (just- accepted):1–21, 2025

    work page 2025

  7. [7]

    Chigansky, P

    P. Chigansky, P. Jagers, and F. C. Klebaner. What can be observed in real time PCR and when does it show? Journal of Mathematical Biology , 76(3):679–695, 2018

    work page 2018

  8. [8]

    Derfel, P

    G. Derfel, P. J. Grabner, and F. Vogl. Asymptotics of the Poincaré functions. In D. Dawson, V. Jakšić, and B. Vainberg, editors, Probability and Mathematical Physics: A Volume in Honor of Stanislav Molchanov , volume 42 of CRM Proceedings and Lecture Notes , pages 113–130. Centre de Recherches Mathématiques, Montréal, 2007

    work page 2007

  9. [9]

    M. S. Dubuc. La densité de la loi-limite d’un processus en cascade expansif. Z. Wahrscheinlichkeit- stheorie und Verw. Gebiete , 19:281–290, 1971

    work page 1971

  10. [10]

    P. L. Duren. Theory of H p Spaces, volume 38. Academic press, 1970. 24 0 20 40 60 80 -1.5 -1 -0.5 0 0.5 1 1.5 Coeffi cients of ϕ c=0.7 Newton c=0.7 forward c=0.9 Newton c=0.9 forward Figure 9: Approximations of the coefficients of φ obtained by Newton’s method and the forward substi- tution method for the two linear fractional examples considered in Section ...

    work page 1970

  11. [11]

    Edelman and G

    A. Edelman and G. Strang. Pascal matrices. The American Mathematical Monthly , 111(3):189–197, 2004

    work page 2004

  12. [12]

    Fernley and E

    J. Fernley and E. Jacob. A universal right tail upper bound for supercritical Galton-Watson processes with bounded offspring. Statist. Probab. Lett. , 209:Paper No. 110082, 5, 2024

    work page 2024

  13. [13]

    T. E. Harris. The theory of branching processes , volume Band 119 of Die Grundlehren der mathe- matischen Wissenschaften . Springer-Verlag, Berlin; Prentice Hall, Inc., Englewood Cliffs, NJ, 1963

    work page 1963

  14. [14]

    Hautphenne and S

    S. Hautphenne and S. Massei. A low-rank technique for computing the quasi-stationary distribu- tion of subcritical Galton–Watson processes. SIAM Journal on Matrix Analysis and Applications , 41(1):29–57, 2020

    work page 2020

  15. [15]

    N. J. Higham. Accuracy and stability of numerical algorithms . SIAM, 2002

    work page 2002

  16. [16]

    Kinoshita and B

    Y. Kinoshita and B. Li. Power series composition in near-linear time. In 2024 IEEE 65th Annual Symposium on Foundations of Computer Science (FOCS) , pages 2180–2185. IEEE, 2024

    work page 2024

  17. [17]

    Morris, J

    D. Morris, J. Maclean, and A. J. Black. Computation of random time-shift distributions for stochastic population models. Journal of Mathematical Biology , 89(3):33, 2024

    work page 2024

  18. [18]

    Poincaré

    H. Poincaré. Sur une classe nouvelle de transcendantes uniformes. Journal de Mathématiques pures et appliquées , 6:313–365, 1890

    work page

  19. [19]

    S. B. Provost and M. Jiang. Orthogonal polynomial density estimates: alternative representation and degree selection. Int. J. Comput. Math. Sci. , 6:17–24, 2012

    work page 2012

  20. [20]

    G. Valiron. Fonctions analytiques. Presses Universitaires de France, Paris, 1954. 25

    work page 1954