← All changes | libraries/vendor/PhpSpreadsheet/Calculation/Statistical/Distributions/GammaBase.php +119 -11 3.3.1 → 3.4 View file →

		@@ -10,9 +10,9 @@
10	10	private const LOG_GAMMA_X_MAX_VALUE = 2.55e305;
11	11
12	12	private const EPS = 2.22e-16;
13	13
14		- private const MAX_VALUE = 1.2e308;
	14	+ protected const MAX_VALUE = 1.2e308;
15	15
16	16	private const SQRT2PI = 2.5066282746310005024157652848110452530069867406099;
17	17
18	18	private const MAX_ITERATIONS = 256;
		@@ -19,12 +19,21 @@
19	19
20	20	protected static function calculateDistribution(float $value, float $a, float $b, bool $cumulative): float
21	21	{
22	22	if ($cumulative) {
23		- return self::incompleteGamma($a, $value / $b) / self::gammaValue($a);
	23	+ return self::regularizedGammaP($a, $value / $b);
24	24	}
25	25
26		- return (1 / ($b ** $a * self::gammaValue($a))) * $value ** ($a - 1) * exp(0 - ($value / $b));
	26	+ if ($value == 0.0) {
	27	+ if ($a == 1.0) {
	28	+ return 1.0 / $b;
	29	+ }
	30	+
	31	+ return ($a < 1.0) ? INF : 0.0;
	32	+ }
	33	+
	34	+ // Log domain, so large shape parameters cannot overflow Gamma(a).
	35	+ return exp(($a - 1.0) * log($value) - $value / $b - $a * log($b) - self::logGamma($a));
27	36	}
28	37
29	38	/** @return float|string */
30	39	protected static function calculateInverse(float $probability, float $alpha, float $beta)
		@@ -31,8 +40,26 @@
31	40	{
32	41	$xLo = 0;
33	42	$xHi = $alpha * $beta * 5;
34	43
	44	+ // Extend the upper bound while it does not yet bracket the root, so a
	45	+ // tail quantile beyond alpha*beta*5 is no longer clamped to it. Stop if
	46	+ // the CDF stops increasing (series approximation past its usable range)
	47	+ // and keep the original bound rather than expanding into that region.
	48	+ $xHiBase = $xHi;
	49	+ $cdfHi = self::calculateDistribution($xHi, $alpha, $beta, true);
	50	+ while ($cdfHi < $probability) {
	51	+ $xHiNext = $xHi * 2;
	52	+ $cdfNext = self::calculateDistribution($xHiNext, $alpha, $beta, true);
	53	+ if ($cdfNext <= $cdfHi) {
	54	+ $xHi = $xHiBase;
	55	+
	56	+ break;
	57	+ }
	58	+ $xHi = $xHiNext;
	59	+ $cdfHi = $cdfNext;
	60	+ }
	61	+
35	62	$dx = 1024;
36	63	$x = $xNew = 1;
37	64	$i = 0;
38	65
		@@ -77,19 +104,100 @@
77	104	// Implementation of the incomplete Gamma function
78	105	//
79	106	public static function incompleteGamma(float $a, float $x): float
80	107	{
81		- static $max = 32;
82		- $summer = 0;
83		- for ($n = 0; $n <= $max; ++$n) {
84		- $divisor = $a;
85		- for ($i = 1; $i <= $n; ++$i) {
86		- $divisor *= ($a + $i);
	108	+ // Unregularized lower incomplete gamma; kept for backward compatibility.
	109	+ return self::regularizedGammaP($a, $x) * self::gammaValue($a);
	110	+ }
	111	+
	112	+ /**
	113	+ * Regularized lower incomplete gamma P(a,x) = gamma(a,x) / Gamma(a).
	114	+ * Series for x < a+1, else the complement of the continued fraction.
	115	+ */
	116	+ public static function regularizedGammaP(float $a, float $x): float
	117	+ {
	118	+ if ($x <= 0.0 || $a <= 0.0) {
	119	+ return 0.0;
	120	+ }
	121	+ if ($x < $a + 1.0) {
	122	+ return self::gammaSeries($a, $x);
	123	+ }
	124	+
	125	+ return 1.0 - self::gammaContinuedFraction($a, $x);
	126	+ }
	127	+
	128	+ /**
	129	+ * Regularized upper incomplete gamma Q(a,x) = 1 - P(a,x).
	130	+ * Continued fraction for x >= a+1 keeps the right tail free of cancellation.
	131	+ */
	132	+ public static function regularizedGammaQ(float $a, float $x): float
	133	+ {
	134	+ if ($x <= 0.0 || $a <= 0.0) {
	135	+ return 1.0;
	136	+ }
	137	+ if ($x < $a + 1.0) {
	138	+ return 1.0 - self::gammaSeries($a, $x);
	139	+ }
	140	+
	141	+ return self::gammaContinuedFraction($a, $x);
	142	+ }
	143	+
	144	+ // Near x ~ a both expansions need O(sqrt(a)) terms to reach EPS.
	145	+ private static function incompleteGammaIterations(float $a): int
	146	+ {
	147	+ return max(self::MAX_ITERATIONS, (int) ceil(10.0 * sqrt($a)));
	148	+ }
	149	+
	150	+ // P(a,x) by its series representation (Numerical Recipes gser).
	151	+ private static function gammaSeries(float $a, float $x): float
	152	+ {
	153	+ $maxIterations = self::incompleteGammaIterations($a);
	154	+ $gln = self::logGamma($a);
	155	+ $ap = $a;
	156	+ $sum = 1.0 / $a;
	157	+ $del = $sum;
	158	+ for ($i = 1; $i <= $maxIterations; ++$i) {
	159	+ ++$ap;
	160	+ $del *= $x / $ap;
	161	+ $sum += $del;
	162	+ if (abs($del) < abs($sum) * self::EPS) {
	163	+ break;
87	164	}
88		- $summer += ($x ** $n / $divisor);
89	165	}
90	166
91		- return $x ** $a * exp(0 - $x) * $summer;
	167	+ return $sum * exp(-$x + $a * log($x) - $gln);
	168	+ }
	169	+
	170	+ // Q(a,x) by its continued fraction representation (Numerical Recipes gcf).
	171	+ private static function gammaContinuedFraction(float $a, float $x): float
	172	+ {
	173	+ $maxIterations = self::incompleteGammaIterations($a);
	174	+ $fpMin = 1.0e-300;
	175	+ $gln = self::logGamma($a);
	176	+ $b = $x + 1.0 - $a;
	177	+ $c = 1.0 / $fpMin;
	178	+ $d = 1.0 / $b;
	179	+ $h = $d;
	180	+ for ($i = 1; $i <= $maxIterations; ++$i) {
	181	+ $an = -$i * ($i - $a);
	182	+ $b += 2.0;
	183	+ $d = $an * $d + $b;
	184	+ if (abs($d) < $fpMin) {
	185	+ $d = $fpMin;
	186	+ }
	187	+ $c = $b + $an / $c;
	188	+ if (abs($c) < $fpMin) {
	189	+ $c = $fpMin;
	190	+ }
	191	+ $d = 1.0 / $d;
	192	+ $del = $d * $c;
	193	+ $h *= $del;
	194	+ if (abs($del - 1.0) < self::EPS) {
	195	+ break;
	196	+ }
	197	+ }
	198	+
	199	+ return $h * exp(-$x + $a * log($x) - $gln);
92	200	}
93	201
94	202	private const GAMMA_VALUE_P0 = 1.000000000190015;
95	203	private const GAMMA_VALUE_P = [