Convergence Methods

Computing moments over infinite discrete support or unbounded continuous intervals requires more than naive summation or integration. ToFUL applies a layered convergence strategy that selects the most effective method for each situation automatically.

The Problem 

For an infinite discrete series, we cannot sum infinitely many terms. We must truncate at some \(N\) and estimate the tail. A naive approach — summing exactly \(N\) terms and stopping — is fragile:

It wastes computation on slowly converging series.
It has no way to detect whether more terms are needed.
It gives no error estimate.

Convergence acceleration algorithms extract accurate limit estimates from the partial sums of a sequence using far fewer terms than naive truncation would require.

Discrete Series: The Convergence Cascade 

ToFUL evaluates all PMF terms over the support using vectorised NumPy, forming partial sums \(S_N = \sum_{k=0}^{N} a_k\). It then applies the following cascade, stopping at the first method that succeeds.

1. Term Magnitude Test 

If the absolute value of each of the last 10 terms is below the convergence tolerance (default \(10^{-12}\)), the partial sum is accepted as converged.

This is the fastest check and catches most well-behaved series (Geometric, Poisson, Binomial with finite support).

2. Wynn ε-Algorithm 

The Wynn epsilon algorithm applies to the sequence of partial sums \(S_0, S_1, \ldots, S_N\). It builds a two-dimensional table via:

\[\varepsilon^{(0)}_n = S_n, \qquad \varepsilon^{(-1)}_n = 0\]

\[\varepsilon^{(k+1)}_n = \varepsilon^{(k-1)}_{n+1} + \frac{1}{\varepsilon^{(k)}_{n+1} - \varepsilon^{(k)}_n}\]

The even-column diagonal entries \(\varepsilon^{(2k)}_0\) converge to the true limit. This method is highly effective for algebraically converging series and is used internally by computer algebra systems such as Mathematica.

The result is accepted when two consecutive even-column estimates agree to within the tolerance and the partial sum matches the Wynn estimate.

3. Aitken Δ² Method 

A lightweight extrapolation using three consecutive partial sums:

\[S^* = S_{n-2} - \frac{(S_{n-1} - S_{n-2})^2}{S_n - 2S_{n-1} + S_{n-2}}\]

The result is accepted when \(|S_N - S^*|\) is below the tolerance. This is the Aitken Δ² method, also called the Shanks transformation of order 1.

4. Cohen-Villegas-Zagier Algorithm 

For alternating series (terms alternate in sign and decrease in absolute value), ToFUL applies the Cohen-Villegas-Zagier algorithm. It computes optimal Chebyshev weights \(d_k\) and evaluates:

\[S = \frac{1}{d_n} \sum_{k=0}^{n-1} (-1)^k (d_k - d_n)\, f(k)\]

This achieves near-machine-precision accuracy in exactly \(n\) term evaluations — far fewer than other methods for alternating series.

5. Geometric Ratio Bound 

If the ratio of consecutive absolute terms converges to \(r < 1\), the tail from term \(N\) onwards is bounded by a geometric series:

\[\text{tail} \approx a_N \cdot \frac{r}{1 - r}\]

This estimate is added to the partial sum. If the ratio is \(r > 1\), the series is flagged as divergent.

6. Partial Sum Fallback 

If none of the above methods achieve convergence, the raw partial sum is returned with a “convergence uncertain” flag. The Convergence tab shows this as a red badge. Increase Max series terms or Calc precision to improve this.

Continuous Integration: Quadrature Selection 

For continuous distributions, ToFUL selects the integration method based on the domain.

Gauss-Laguerre Quadrature 

Domain: \([0, +\infty)\)

The \(n\)-point Gauss-Laguerre rule uses nodes \(x_i\) and weights \(w_i\) satisfying:

\[\int_0^\infty f(x)\, e^{-x}\, dx \approx \sum_{i=1}^n w_i\, f(x_i)\]

To integrate \(f(x)\) (without the \(e^{-x}\) weight), we evaluate \(f(x_i) \cdot e^{x_i}\) at the nodes:

\[\int_0^\infty f(x)\, dx \approx \sum_{i=1}^n w_i\, f(x_i)\, e^{x_i}\]

With \(n = 64\) nodes, this is exact for polynomials of degree up to \(2n - 1 = 127\) and highly accurate for smooth exponentially decaying functions. The result is cross-checked against SciPy adaptive quad; the closer of the two estimates is returned.

Gauss-Hermite Quadrature 

Domain: \((-\infty, +\infty)\)

The \(n\)-point Gauss-Hermite rule targets integrals weighted by \(e^{-x^2}\). For unweighted integrals:

\[\int_{-\infty}^{\infty} f(x)\, dx \approx \sum_{i=1}^n w_i\, f(x_i)\, e^{x_i^2}\]

With \(n = 64\) nodes. Cross-checked against adaptive quad.

Adaptive Gauss-Kronrod (SciPy `quad`)

Domain: Any finite interval, or semi-infinite domains not starting at 0.

SciPy’s quad uses adaptive Gauss-Kronrod rules (15-point / 7-point pairs) with automatic subdivision. Tolerances are set to \(10^{-(\text{prec}+4)}\) where prec is the Calc precision setting.

mpmath Tanh-Sinh Quadrature 

Triggered by: Calc precision ≥ 13

The tanh-sinh (double-exponential) method applies the substitution \(x = \tanh(\pi \sinh t / 2)\), which concentrates quadrature nodes near endpoint singularities. It achieves exponential convergence for a wide class of integrands, including those with endpoint singularities, and supports arbitrary precision via mpmath.

SymPy Symbolic Integration 

Triggered by: Expressions without if/else conditionals

Before any numerical integration, ToFUL attempts to compute the integral symbolically via SymPy. If successful, the result is exact (no numerical error) and the method shown is sympy-exact. This fires most reliably for classical distributions like Normal, Gamma, and Beta.

Richardson Error Estimation 

For the specialised Gauss rules, ToFUL estimates the integration error by comparing the \(n = 64\) node result against a \(n = 32\) node result (Richardson extrapolation). The difference is reported as the error estimate in the Convergence tab info string.