| declare(strict_types=1); | |
| namespace Dochia; | |
| /** | |
| * Normal/lognormal numeric helpers (no extensions, no dependencies). | |
| * | |
| * erf uses Abramowitz & Stegun 7.1.26 (|error| < 1.5e-7), which is ample for | |
| * day-level probabilities and quantiles; golden tests compare with the Python | |
| * reference at a tolerant epsilon. | |
| */ | |
| final class Erf | |
| { | |
| public static function erf(float $x): float | |
| { | |
| $sign = $x < 0 ? -1.0 : 1.0; | |
| $x = abs($x); | |
| $a1 = 0.254829592; | |
| $a2 = -0.284496736; | |
| $a3 = 1.421413741; | |
| $a4 = -1.453152027; | |
| $a5 = 1.061405429; | |
| $p = 0.3275911; | |
| $t = 1.0 / (1.0 + $p * $x); | |
| $y = 1.0 - ((((($a5 * $t + $a4) * $t + $a3) * $t + $a2) * $t + $a1) * $t) * exp(-$x * $x); | |
| return $sign * $y; | |
| } | |
| public static function normCdf(float $z): float | |
| { | |
| return 0.5 * (1.0 + self::erf($z / M_SQRT2)); | |
| } | |
| public static function normPdf(float $z): float | |
| { | |
| return exp(-0.5 * $z * $z) / sqrt(2.0 * M_PI); | |
| } | |
| /** Inverse normal CDF by bisection (matches the Python reference to ~1e-12). */ | |
| public static function normPpf(float $p): float | |
| { | |
| if ($p <= 0.0) { | |
| return -INF; | |
| } | |
| if ($p >= 1.0) { | |
| return INF; | |
| } | |
| $lo = -12.0; | |
| $hi = 12.0; | |
| for ($i = 0; $i < 200; $i++) { | |
| $mid = 0.5 * ($lo + $hi); | |
| if (self::normCdf($mid) < $p) { | |
| $lo = $mid; | |
| } else { | |
| $hi = $mid; | |
| } | |
| } | |
| return 0.5 * ($lo + $hi); | |
| } | |
| } | |