diff --git a/src/Functions/Special.php b/src/Functions/Special.php index d2e8845c5..b539b9936 100644 --- a/src/Functions/Special.php +++ b/src/Functions/Special.php @@ -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 * @@ -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; @@ -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; } /** @@ -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) / $√π; } /** @@ -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); } /** diff --git a/tests/Functions/Special/ErrorFunctionTest.php b/tests/Functions/Special/ErrorFunctionTest.php index 4e226366d..a080a5877 100644 --- a/tests/Functions/Special/ErrorFunctionTest.php +++ b/tests/Functions/Special/ErrorFunctionTest.php @@ -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"); + } + } } diff --git a/tests/Probability/Distribution/Continuous/NormalTest.php b/tests/Probability/Distribution/Continuous/NormalTest.php index 985230495..a99b06ff6 100644 --- a/tests/Probability/Distribution/Continuous/NormalTest.php +++ b/tests/Probability/Distribution/Continuous/NormalTest.php @@ -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], + ]; + } }