tablepress / libraries / vendor / PhpSpreadsheet / Calculation / Statistical / Distributions / GammaBase.php

GammaBase.php in TablePress – Tables in WordPress made easy 3.4, at libraries/vendor/PhpSpreadsheet/Calculation/Statistical/Distributions/GammaBase.php

491 lines 12.3 KB

Raw Download Zip

1	<?php
2
3	namespace TablePress\PhpOffice\PhpSpreadsheet\Calculation\Statistical\Distributions;
4
5	use TablePress\PhpOffice\PhpSpreadsheet\Calculation\Functions;
6	use TablePress\PhpOffice\PhpSpreadsheet\Calculation\Information\ExcelError;
7
8	abstract class GammaBase
9	{
10	private const LOG_GAMMA_X_MAX_VALUE = 2.55e305;
11
12	private const EPS = 2.22e-16;
13
14	protected const MAX_VALUE = 1.2e308;
15
16	private const SQRT2PI = 2.5066282746310005024157652848110452530069867406099;
17
18	private const MAX_ITERATIONS = 256;
19
20	protected static function calculateDistribution(float $value, float $a, float $b, bool $cumulative): float
21	{
22	if ($cumulative) {
23	return self::regularizedGammaP($a, $value / $b);
24	}
25
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));
36	}
37
38	/** @return float\|string */
39	protected static function calculateInverse(float $probability, float $alpha, float $beta)
40	{
41	$xLo = 0;
42	$xHi = $alpha * $beta * 5;
43
44	// Extend the upper bound while it does not yet bracket the root, so a
45	// tail quantile beyond alphabeta5 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
62	$dx = 1024;
63	$x = $xNew = 1;
64	$i = 0;
65
66	while ((abs($dx) > Functions::PRECISION) && (++$i <= self::MAX_ITERATIONS)) {
67	// Apply Newton-Raphson step
68	$result = self::calculateDistribution($x, $alpha, $beta, true);
69	$error = $result - $probability;
70
71	if ($error == 0.0) {
72	$dx = 0;
73	} elseif ($error < 0.0) {
74	$xLo = $x;
75	} else {
76	$xHi = $x;
77	}
78
79	$pdf = self::calculateDistribution($x, $alpha, $beta, false);
80	// Avoid division by zero
81	if ($pdf !== 0.0) {
82	$dx = $error / $pdf;
83	$xNew = $x - $dx;
84	}
85
86	// If the NR fails to converge (which for example may be the
87	// case if the initial guess is too rough) we apply a bisection
88	// step to determine a more narrow interval around the root.
89	if (($xNew < $xLo) \|\| ($xNew > $xHi) \|\| ($pdf == 0.0)) {
90	$xNew = ($xLo + $xHi) / 2;
91	$dx = $xNew - $x;
92	}
93	$x = $xNew;
94	}
95
96	if ($i === self::MAX_ITERATIONS) {
97	return ExcelError::NA();
98	}
99
100	return $x;
101	}
102
103	//
104	// Implementation of the incomplete Gamma function
105	//
106	public static function incompleteGamma(float $a, float $x): float
107	{
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;
164	}
165	}
166
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);
200	}
201
202	private const GAMMA_VALUE_P0 = 1.000000000190015;
203	private const GAMMA_VALUE_P = [
204	1 => 76.18009172947146,
205	2 => -86.50532032941677,
206	3 => 24.01409824083091,
207	4 => -1.231739572450155,
208	5 => 1.208650973866179e-3,
209	6 => -5.395239384953e-6,
210	];
211
212	//
213	// Implementation of the Gamma function
214	//
215	public static function gammaValue(float $value): float
216	{
217	if ($value == 0.0) {
218	return 0;
219	}
220
221	$y = $x = $value;
222	$tmp = $x + 5.5;
223	$tmp -= ($x + 0.5) * log($tmp);
224
225	$summer = self::GAMMA_VALUE_P0;
226	for ($j = 1; $j <= 6; ++$j) {
227	$summer += (self::GAMMA_VALUE_P[$j] / ++$y);
228	}
229
230	return exp(0 - $tmp + log(self::SQRT2PI * $summer / $x));
231	}
232
233	private const LG_D1 = -0.5772156649015328605195174;
234
235	private const LG_D2 = 0.4227843350984671393993777;
236
237	private const LG_D4 = 1.791759469228055000094023;
238
239	private const LG_P1 = [
240	4.945235359296727046734888,
241	201.8112620856775083915565,
242	2290.838373831346393026739,
243	11319.67205903380828685045,
244	28557.24635671635335736389,
245	38484.96228443793359990269,
246	26377.48787624195437963534,
247	7225.813979700288197698961,
248	];
249
250	private const LG_P2 = [
251	4.974607845568932035012064,
252	542.4138599891070494101986,
253	15506.93864978364947665077,
254	184793.2904445632425417223,
255	1088204.76946882876749847,
256	3338152.967987029735917223,
257	5106661.678927352456275255,
258	3074109.054850539556250927,
259	];
260
261	private const LG_P4 = [
262	14745.02166059939948905062,
263	2426813.369486704502836312,
264	121475557.4045093227939592,
265	2663432449.630976949898078,
266	29403789566.34553899906876,
267	170266573776.5398868392998,
268	492612579337.743088758812,
269	560625185622.3951465078242,
270	];
271
272	private const LG_Q1 = [
273	67.48212550303777196073036,
274	1113.332393857199323513008,
275	7738.757056935398733233834,
276	27639.87074403340708898585,
277	54993.10206226157329794414,
278	61611.22180066002127833352,
279	36351.27591501940507276287,
280	8785.536302431013170870835,
281	];
282
283	private const LG_Q2 = [
284	183.0328399370592604055942,
285	7765.049321445005871323047,
286	133190.3827966074194402448,
287	1136705.821321969608938755,
288	5267964.117437946917577538,
289	13467014.54311101692290052,
290	17827365.30353274213975932,
291	9533095.591844353613395747,
292	];
293
294	private const LG_Q4 = [
295	2690.530175870899333379843,
296	639388.5654300092398984238,
297	41355999.30241388052042842,
298	1120872109.61614794137657,
299	14886137286.78813811542398,
300	101680358627.2438228077304,
301	341747634550.7377132798597,
302	446315818741.9713286462081,
303	];
304
305	private const LG_C = [
306	-0.001910444077728,
307	8.4171387781295e-4,
308	-5.952379913043012e-4,
309	7.93650793500350248e-4,
310	-0.002777777777777681622553,
311	0.08333333333333333331554247,
312	0.0057083835261,
313	];
314
315	// Rough estimate of the fourth root of logGamma_xBig
316	private const LG_FRTBIG = 2.25e76;
317
318	private const PNT68 = 0.6796875;
319
320	// Function cache for logGamma
321
322	private static float $logGammaCacheResult = 0.0;
323
324	private static float $logGammaCacheX = 0.0;
325
326	/**
327	* logGamma function.
328	*
329	* Original author was Jaco van Kooten. Ported to PHP by Paul Meagher.
330	*
331	* The natural logarithm of the gamma function. <br />
332	* Based on public domain NETLIB (Fortran) code by W. J. Cody and L. Stoltz <br />
333	* Applied Mathematics Division <br />
334	* Argonne National Laboratory <br />
335	* Argonne, IL 60439 <br />
336	* <p>
337	* References:
338	* <ol>
339	* <li>W. J. Cody and K. E. Hillstrom, 'Chebyshev Approximations for the Natural
340	* Logarithm of the Gamma Function,' Math. Comp. 21, 1967, pp. 198-203.</li>
341	* <li>K. E. Hillstrom, ANL/AMD Program ANLC366S, DGAMMA/DLGAMA, May, 1969.</li>
342	* <li>Hart, Et. Al., Computer Approximations, Wiley and sons, New York, 1968.</li>
343	* </ol>
344	* </p>
345	* <p>
346	* From the original documentation:
347	* </p>
348	* <p>
349	* This routine calculates the LOG(GAMMA) function for a positive real argument X.
350	* Computation is based on an algorithm outlined in references 1 and 2.
351	* The program uses rational functions that theoretically approximate LOG(GAMMA)
352	* to at least 18 significant decimal digits. The approximation for X > 12 is from
353	* reference 3, while approximations for X < 12.0 are similar to those in reference
354	* 1, but are unpublished. The accuracy achieved depends on the arithmetic system,
355	* the compiler, the intrinsic functions, and proper selection of the
356	* machine-dependent constants.
357	* </p>
358	* <p>
359	* Error returns: <br />
360	* The program returns the value XINF for X .LE. 0.0 or when overflow would occur.
361	* The computation is believed to be free of underflow and overflow.
362	* </p>
363	*
364	* @version 1.1
365	*
366	* @author Jaco van Kooten
367	*
368	* @return float MAX_VALUE for x < 0.0 or when overflow would occur, i.e. x > 2.55E305
369	*/
370	public static function logGamma(float $x): float
371	{
372	if ($x == self::$logGammaCacheX) {
373	return self::$logGammaCacheResult;
374	}
375
376	$y = $x;
377	if ($y > 0.0 && $y <= self::LOG_GAMMA_X_MAX_VALUE) {
378	if ($y <= self::EPS) {
379	$res = -log($y);
380	} elseif ($y <= 1.5) {
381	$res = self::logGamma1($y);
382	} elseif ($y <= 4.0) {
383	$res = self::logGamma2($y);
384	} elseif ($y <= 12.0) {
385	$res = self::logGamma3($y);
386	} else {
387	$res = self::logGamma4($y);
388	}
389	} else {
390	// --------------------------
391	// Return for bad arguments
392	// --------------------------
393	$res = self::MAX_VALUE;
394	}
395
396	// ------------------------------
397	// Final adjustments and return
398	// ------------------------------
399	self::$logGammaCacheX = $x;
400	self::$logGammaCacheResult = $res;
401
402	return $res;
403	}
404
405	private static function logGamma1(float $y): float
406	{
407	// ---------------------
408	// EPS .LT. X .LE. 1.5
409	// ---------------------
410	if ($y < self::PNT68) {
411	$corr = -log($y);
412	$xm1 = $y;
413	} else {
414	$corr = 0.0;
415	$xm1 = $y - 1.0;
416	}
417
418	$xden = 1.0;
419	$xnum = 0.0;
420	if ($y <= 0.5 \|\| $y >= self::PNT68) {
421	for ($i = 0; $i < 8; ++$i) {
422	$xnum = $xnum * $xm1 + self::LG_P1[$i];
423	$xden = $xden * $xm1 + self::LG_Q1[$i];
424	}
425
426	return $corr + $xm1 * (self::LG_D1 + $xm1 * ($xnum / $xden));
427	}
428
429	$xm2 = $y - 1.0;
430	for ($i = 0; $i < 8; ++$i) {
431	$xnum = $xnum * $xm2 + self::LG_P2[$i];
432	$xden = $xden * $xm2 + self::LG_Q2[$i];
433	}
434
435	return $corr + $xm2 * (self::LG_D2 + $xm2 * ($xnum / $xden));
436	}
437
438	private static function logGamma2(float $y): float
439	{
440	// ---------------------
441	// 1.5 .LT. X .LE. 4.0
442	// ---------------------
443	$xm2 = $y - 2.0;
444	$xden = 1.0;
445	$xnum = 0.0;
446	for ($i = 0; $i < 8; ++$i) {
447	$xnum = $xnum * $xm2 + self::LG_P2[$i];
448	$xden = $xden * $xm2 + self::LG_Q2[$i];
449	}
450
451	return $xm2 * (self::LG_D2 + $xm2 * ($xnum / $xden));
452	}
453
454	protected static function logGamma3(float $y): float
455	{
456	// ----------------------
457	// 4.0 .LT. X .LE. 12.0
458	// ----------------------
459	$xm4 = $y - 4.0;
460	$xden = -1.0;
461	$xnum = 0.0;
462	for ($i = 0; $i < 8; ++$i) {
463	$xnum = $xnum * $xm4 + self::LG_P4[$i];
464	$xden = $xden * $xm4 + self::LG_Q4[$i];
465	}
466
467	return self::LG_D4 + $xm4 * ($xnum / $xden);
468	}
469
470	protected static function logGamma4(float $y): float
471	{
472	// ---------------------------------
473	// Evaluate for argument .GE. 12.0
474	// ---------------------------------
475	$res = 0.0;
476	if ($y <= self::LG_FRTBIG) {
477	$res = self::LG_C[6];
478	$ysq = $y * $y;
479	for ($i = 0; $i < 6; ++$i) {
480	$res = $res / $ysq + self::LG_C[$i];
481	}
482	$res /= $y;
483	$corr = log($y);
484	$res = $res + log(self::SQRT2PI) - 0.5 * $corr;
485	$res += $y * ($corr - 1.0);
486	}
487
488	return $res;
489	}
490	}
491