Skip to content

Fix erf/erfc catastrophic cancellation producing impossible values (erf > 1, erfc < 0, negative Normal survival) #491

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
gaoflow wants to merge 1 commit into
base: develop
Choose a base branch
from
Open

Fix erf/erfc catastrophic cancellation producing impossible values (erf > 1, erfc < 0, negative Normal survival) #491

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
92 changes: 27 additions & 65 deletions src/Functions/Special.php
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -915,18 +915,16 @@ public static function sigmoid(float $t): float
* erf(x) = ── ⎮ e^(-t²) dt
* √π ⌡₀
*
* Improved implementation with domain-specific algorithms:
* - Small arguments (|x| ≤ 0.01): Taylor series (8 terms, optimized)
* - Medium arguments (0.01 < |x| ≤ 4): Taylor series with convergence (up to 50 terms)
* - Large arguments (|x| > 4): Asymptotic expansion (4 terms)
*
* This implementation prioritizes accuracy (< 1e-12 error) over the classical
* Abramowitz & Stegun 7.1.26 approximation (max error: 1.5e-7).
* Algorithm:
* - Small arguments (|x| ≤ 0.01): 8-term Maclaurin series (avoids 1 - erfc cancellation near 0)
* - Otherwise: erf(x) = sign(x) · (1 - erfc(|x|)), with erfc from the regularized upper
* incomplete gamma identity erfc(x) = Γ(½, x²)/√π
*
* Taylor series: erf(x) = (2/√π) * Σ((-1)^n * x^(2n+1) / (n! * (2n+1)))
* Asymptotic: erf(x) = 1 - erfc(x) where erfc(x) ≈ (e^(-x²) / (x√π)) * (1 - 1/(2x²) + ...)
* The alternating Maclaurin series for erf loses ~11 digits to cancellation near x = 4
* (terms grow to ≈ e^{x²} before decaying), so it is only used for tiny x. Everywhere else
* erf is derived from the numerically stable gamma continued fraction already in this file.
*
* Precision: Better than 1e-12 for all x
* Maclaurin series: erf(x) = (2/√π) * (x - x³/3 + x⁵/10 - x⁷/42 + ...)
*
* @param float $x
*
Expand All @@ -942,7 +940,7 @@ public static function errorFunction(float $x): float

$two／√π = 1.1283791670955125738961589031215451716881;

// Small arguments: Use Taylor series
// Small arguments: Maclaurin series (stable for tiny x, where 1 - erfc would cancel)
// erf(x) = (2/√π) * (x - x³/3 + x⁵/10 - x⁷/42 + x⁹/216 - ...)
if ($ax <= 0.01) {
$x² = $x * $x;
Expand All @@ -959,37 +957,10 @@ public static function errorFunction(float $x): float
return $sum * $two／√π;
}

// Large arguments: For x > 4, compute using erfc for consistency
// erf(x) = 1 - erfc(x), but computed via asymptotic expansion to avoid cancellation
if ($ax > 4.0) {
// For positive x: erf(x) = 1 - erfc(x)
if ($x > 0) {
return 1.0 - self::erfcAsymptoticSeries($x);
} else {
// For negative x: erf(-x) = -erf(x)
return -(1.0 - self::erfcAsymptoticSeries($ax));
}
}

// Medium arguments: Use Taylor series
// erf(x) = (2/√π) * Σ((-1)^n * x^(2n+1) / (n! * (2n+1)))
$x² = $ax * $ax;
$sum = $ax;
$term = $ax;

for ($n = 1; $n <= 50; $n++) {
$term *= -$x² / $n;
$sum += $term / (2 * $n + 1);

// Early exit if converged
if (\abs($term / (2 * $n + 1)) < 1e-15) {
break;
}
}

$result = $sum * $two／√π;
// erf(x) = sign(x) · (1 - erfc(|x|)); erfc via stable Γ(½, x²)/√π
$erf = 1.0 - self::erfcViaIncompleteGamma($ax);

return ($x > 0) ? $result : -$result;
return ($x > 0) ? $erf : -$erf;
}

/**
Expand All @@ -1006,27 +977,21 @@ public static function erf(float $x): float
}

/**
* Asymptotic series for erfc
* erfc(x) ≈ (e^(-x²) / (x√π)) * (1 - 1/(2x²) + 3/(4x⁴) - ...)
* erfc for non-negative x via the regularized upper incomplete gamma function
* erfc(x) = Γ(½, x²) / √π (x ≥ 0)
*
* @param float $x Must be positive and >= 4.0
* Reuses the stable continued-fraction / series machinery of upperIncompleteGamma,
* which does not suffer the cancellation of the alternating erf Maclaurin series.
*
* @param float $x Must be >= 0
*
* @return float
*/
private static function erfcAsymptoticSeries(float $x): float
private static function erfcViaIncompleteGamma(float $x): float
{
$x² = $x * $x;
$x⁴ = $x² * $x²;
$x⁶ = $x⁴ * $x²;

$series = 1.0;
$series -= 1.0 / (2.0 * $x²);
$series += 3.0 / (4.0 * $x⁴);
$series -= 15.0 / (8.0 * $x⁶);
$√π = 1.7724538509055160272981674833411451828;

$factor = \exp(-$x²) / ($x * 1.7724538509055160272981674833411451828); // √π

return $factor * $series;
return self::upperIncompleteGamma(0.5, $x * $x) / $√π;
}

/**
Expand All @@ -1052,18 +1017,15 @@ private static function erfcAsymptoticSeries(float $x): float
*/
public static function complementaryErrorFunction($x): float
{
// For large positive x, use asymptotic expansion to avoid 1 - erf(x) cancellation
if ($x >= 4.0) {
return self::erfcAsymptoticSeries($x);
}
$x = (float)$x;

// For large negative x, erfc(-x) = 2 - erfc(x)
if ($x <= -6.0) {
return 2 - self::erfc(-$x);
// erfc(x) = 2 - erfc(-x) for negative x; keeps the reflection exact
if ($x < 0) {
return 2.0 - self::erfcViaIncompleteGamma(-$x);
}

// Otherwise use erf
return 1 - self::errorFunction($x);
// erfc(x) = Γ(½, x²)/√π, computed directly to avoid 1 - erf(x) cancellation
return self::erfcViaIncompleteGamma($x);
}

/**
Expand Down
77 changes: 77 additions & 0 deletions tests/Functions/Special/ErrorFunctionTest.php
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -840,4 +840,81 @@ public function dataProviderForErfcAbramowitzStegun(): array
[2.0, 0.0046777349810472658],
];
}

/**
* @test erf and erfc match arbitrary-precision reference values across every branch,
* and their odd/reflection symmetries hold
* @dataProvider dataProviderForErfErfcReference
* @param float $x
* @param float $erf reference erf(x)
* @param float $erfc reference erfc(x)
*/
public function testErfErfcReference(float $x, float $erf, float $erfc)
{
// When
$erf_pos = Special::erf($x);
$erfc_pos = Special::erfc($x);
$erf_neg = Special::erf(-$x);
$erfc_neg = Special::erfc(-$x);

// Then values match the reference to a tight relative tolerance
$this->assertEqualsWithDelta($erf, $erf_pos, \abs($erf) * 1e-12 + 1e-15);
$this->assertEqualsWithDelta($erfc, $erfc_pos, \abs($erfc) * 1e-12 + 1e-300);

// And the odd (erf(-x) = -erf(x)) and reflection (erfc(-x) = 2 - erfc(x)) identities hold
$this->assertEqualsWithDelta(-$erf, $erf_neg, \abs($erf) * 1e-12 + 1e-15);
$this->assertEqualsWithDelta(2 - $erfc, $erfc_neg, \abs($erfc) * 1e-12 + 1e-15);
}

/**
* Reference erf/erfc values from mpmath (mpmath.mp.dps = 40, cross-checked at dps 35).
* Grid spans every dispatch branch, densest across 3.0 ≤ x ≤ 4.5 where the previous
* alternating Maclaurin series lost ~11 digits to catastrophic cancellation.
* @return array (x, erf, erfc)
*/
public function dataProviderForErfErfcReference(): array
{
return [
[0.005, 0.0056418488200315504, 0.99435815117996845],
[0.01, 0.011283415555849617, 0.98871658444415038],
[0.05, 0.056371977797016627, 0.94362802220298337],
[0.1, 0.11246291601828490, 0.88753708398171510],
[0.25, 0.27632639016823693, 0.72367360983176307],
[0.5, 0.52049987781304654, 0.47950012218695346],
[0.75, 0.71115563365351513, 0.28884436634648487],
[1.0, 0.84270079294971487, 0.15729920705028513],
[1.5, 0.96610514647531073, 0.033894853524689273],
[2.0, 0.99532226501895273, 0.0046777349810472658],
[2.5, 0.99959304798255504, 0.00040695201744495894],
[3.0, 0.99997790950300141, 2.2090496998585441e-5],
[3.3, 0.99999694229020356, 3.0577097964381652e-6],
[3.5, 0.99999925690162766, 7.4309837234141275e-7],
[3.7, 0.99999983284894209, 1.6715105790914598e-7],
[3.9, 0.99999996520775140, 3.4792248597231767e-8],
[3.99, 0.99999998326078864, 1.6739211364520814e-8],
[4.0, 0.99999998458274210, 1.5417257900280019e-8],
[4.1, 0.99999999329997235, 6.7000276540849184e-9],
[4.5, 0.99999999980338396, 1.9661604415428875e-10],
[5.0, 0.99999999999846254, 1.5374597944280349e-12],
[5.5, 0.99999999999999264, 7.3578479179743981e-15],
[6.0, 0.99999999999999998, 2.1519736712498913e-17],
[7.0, 1.0000000000000000, 4.1838256077794144e-23],
[8.0, 1.0000000000000000, 1.1224297172982927e-29],
];
}

/**
* @test erf and erfc never leave their mathematical bounds on a dense sweep
* (regression: cancellation near x = 4 produced erf > 1 and erfc < 0)
*/
public function testErfErfcStayWithinBounds()
{
for ($x = 0.0; $x <= 8.0; $x += 0.01) {
// erf(x) ∈ [-1, 1], erfc(x) ∈ [0, 2] for all real x
$this->assertLessThanOrEqual(1.0, Special::erf($x), "erf($x) must not exceed 1");
$this->assertGreaterThanOrEqual(-1.0, Special::erf(-$x), "erf(-$x) must not be below -1");
$this->assertGreaterThanOrEqual(0.0, Special::erfc($x), "erfc($x) must not be negative");
$this->assertLessThanOrEqual(2.0, Special::erfc(-$x), "erfc(-$x) must not exceed 2");
}
}
}
40 changes: 40 additions & 0 deletions tests/Probability/Distribution/Continuous/NormalTest.php
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -662,4 +662,44 @@ public function testRand()
}
}
}

/**
* @test upper-tail survival probability 1 - cdf(z) is never negative
* (regression: erf cancellation drove the standard-normal tail below zero for z ≈ 5-5.7)
* @dataProvider dataProviderForUpperTailSurvival
* @param float $z
* @param float $expectedSurvival reference 1 - Φ(z) from mpmath (dps 40)
*/
public function testUpperTailSurvivalIsNonNegative(float $z, float $expectedSurvival)
{
// Given
$normal = new Normal(0, 1);

// When
$survival = 1 - $normal->cdf($z);

// Then a survival probability can never be negative
$this->assertGreaterThanOrEqual(0.0, $survival, "survival 1 - cdf($z) must be >= 0");

// And it matches the reference magnitude. 1 - cdf(z) subtracts a value near 1, so its own
// rounding caps the deep tail at ~1e-16 absolute; the floor accounts for that, not for erf.
$this->assertEqualsWithDelta($expectedSurvival, $survival, $expectedSurvival * 1e-6 + 2e-15);
}

/**
* Reference standard-normal survival 1 - Φ(z) = erfc(z/√2)/2 from mpmath (dps 40).
* @return array (z, survival)
*/
public function dataProviderForUpperTailSurvival(): array
{
return [
[4.5, 3.3976731247300604e-6],
[5.0, 2.8665157187919391e-7],
[5.6, 1.0717590258310929e-8],
[5.643, 8.355617331197473e-9],
[6.0, 9.8658764503769814e-10],
[6.5, 4.0160005838591178e-11],
[7.0, 1.279812543885835e-12],
];
}
}
Loading