File size: 1,593 Bytes
e215d3c | 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 | <?php
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);
}
}
|