|
1
|
|
|
<?php |
|
2
|
|
|
|
|
3
|
|
|
namespace PhpOffice\PhpSpreadsheet\Calculation\Statistical\Distributions; |
|
4
|
|
|
|
|
5
|
|
|
use PhpOffice\PhpSpreadsheet\Calculation\Functions; |
|
6
|
|
|
use 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
|
|
|
private const MAX_VALUE = 1.2e308; |
|
15
|
|
|
|
|
16
|
|
|
private const SQRT2PI = 2.5066282746310005024157652848110452530069867406099; |
|
17
|
|
|
|
|
18
|
|
|
private const MAX_ITERATIONS = 256; |
|
19
|
|
|
|
|
20
|
11 |
|
protected static function calculateDistribution(float $value, float $a, float $b, bool $cumulative): float |
|
21
|
|
|
{ |
|
22
|
11 |
|
if ($cumulative) { |
|
23
|
8 |
|
return self::incompleteGamma($a, $value / $b) / self::gammaValue($a); |
|
24
|
|
|
} |
|
25
|
|
|
|
|
26
|
7 |
|
return (1 / ($b ** $a * self::gammaValue($a))) * $value ** ($a - 1) * exp(0 - ($value / $b)); |
|
27
|
|
|
} |
|
28
|
|
|
|
|
29
|
|
|
/** @return float|string */ |
|
30
|
4 |
|
protected static function calculateInverse(float $probability, float $alpha, float $beta) |
|
31
|
|
|
{ |
|
32
|
4 |
|
$xLo = 0; |
|
33
|
4 |
|
$xHi = $alpha * $beta * 5; |
|
34
|
|
|
|
|
35
|
4 |
|
$dx = 1024; |
|
36
|
4 |
|
$x = $xNew = 1; |
|
37
|
4 |
|
$i = 0; |
|
38
|
|
|
|
|
39
|
4 |
|
while ((abs($dx) > Functions::PRECISION) && (++$i <= self::MAX_ITERATIONS)) { |
|
40
|
|
|
// Apply Newton-Raphson step |
|
41
|
4 |
|
$result = self::calculateDistribution($x, $alpha, $beta, true); |
|
42
|
4 |
|
$error = $result - $probability; |
|
43
|
|
|
|
|
44
|
4 |
|
if ($error == 0.0) { |
|
45
|
3 |
|
$dx = 0; |
|
46
|
4 |
|
} elseif ($error < 0.0) { |
|
47
|
4 |
|
$xLo = $x; |
|
48
|
|
|
} else { |
|
49
|
4 |
|
$xHi = $x; |
|
50
|
|
|
} |
|
51
|
|
|
|
|
52
|
4 |
|
$pdf = self::calculateDistribution($x, $alpha, $beta, false); |
|
53
|
|
|
// Avoid division by zero |
|
54
|
4 |
|
if ($pdf !== 0.0) { |
|
55
|
4 |
|
$dx = $error / $pdf; |
|
56
|
4 |
|
$xNew = $x - $dx; |
|
57
|
|
|
} |
|
58
|
|
|
|
|
59
|
|
|
// If the NR fails to converge (which for example may be the |
|
60
|
|
|
// case if the initial guess is too rough) we apply a bisection |
|
61
|
|
|
// step to determine a more narrow interval around the root. |
|
62
|
4 |
|
if (($xNew < $xLo) || ($xNew > $xHi) || ($pdf == 0.0)) { |
|
63
|
4 |
|
$xNew = ($xLo + $xHi) / 2; |
|
64
|
4 |
|
$dx = $xNew - $x; |
|
65
|
|
|
} |
|
66
|
4 |
|
$x = $xNew; |
|
67
|
|
|
} |
|
68
|
|
|
|
|
69
|
4 |
|
if ($i === self::MAX_ITERATIONS) { |
|
70
|
|
|
return ExcelError::NA(); |
|
71
|
|
|
} |
|
72
|
|
|
|
|
73
|
4 |
|
return $x; |
|
74
|
|
|
} |
|
75
|
|
|
|
|
76
|
|
|
// |
|
77
|
|
|
// Implementation of the incomplete Gamma function |
|
78
|
|
|
// |
|
79
|
37 |
|
public static function incompleteGamma(float $a, float $x): float |
|
80
|
|
|
{ |
|
81
|
37 |
|
static $max = 32; |
|
82
|
37 |
|
$summer = 0; |
|
83
|
37 |
|
for ($n = 0; $n <= $max; ++$n) { |
|
84
|
37 |
|
$divisor = $a; |
|
85
|
37 |
|
for ($i = 1; $i <= $n; ++$i) { |
|
86
|
37 |
|
$divisor *= ($a + $i); |
|
87
|
|
|
} |
|
88
|
37 |
|
$summer += ($x ** $n / $divisor); |
|
89
|
|
|
} |
|
90
|
|
|
|
|
91
|
37 |
|
return $x ** $a * exp(0 - $x) * $summer; |
|
92
|
|
|
} |
|
93
|
|
|
|
|
94
|
|
|
private const GAMMA_VALUE_P0 = 1.000000000190015; |
|
95
|
|
|
private const GAMMA_VALUE_P = [ |
|
96
|
|
|
1 => 76.18009172947146, |
|
97
|
|
|
2 => -86.50532032941677, |
|
98
|
|
|
3 => 24.01409824083091, |
|
99
|
|
|
4 => -1.231739572450155, |
|
100
|
|
|
5 => 1.208650973866179e-3, |
|
101
|
|
|
6 => -5.395239384953e-6, |
|
102
|
|
|
]; |
|
103
|
|
|
|
|
104
|
|
|
// |
|
105
|
|
|
// Implementation of the Gamma function |
|
106
|
|
|
// |
|
107
|
78 |
|
public static function gammaValue(float $value): float |
|
108
|
|
|
{ |
|
109
|
78 |
|
if ($value == 0.0) { |
|
110
|
|
|
return 0; |
|
111
|
|
|
} |
|
112
|
|
|
|
|
113
|
78 |
|
$y = $x = $value; |
|
114
|
78 |
|
$tmp = $x + 5.5; |
|
115
|
78 |
|
$tmp -= ($x + 0.5) * log($tmp); |
|
116
|
|
|
|
|
117
|
78 |
|
$summer = self::GAMMA_VALUE_P0; |
|
118
|
78 |
|
for ($j = 1; $j <= 6; ++$j) { |
|
119
|
78 |
|
$summer += (self::GAMMA_VALUE_P[$j] / ++$y); |
|
120
|
|
|
} |
|
121
|
|
|
|
|
122
|
78 |
|
return exp(0 - $tmp + log(self::SQRT2PI * $summer / $x)); |
|
123
|
|
|
} |
|
124
|
|
|
|
|
125
|
|
|
private const LG_D1 = -0.5772156649015328605195174; |
|
126
|
|
|
|
|
127
|
|
|
private const LG_D2 = 0.4227843350984671393993777; |
|
128
|
|
|
|
|
129
|
|
|
private const LG_D4 = 1.791759469228055000094023; |
|
130
|
|
|
|
|
131
|
|
|
private const LG_P1 = [ |
|
132
|
|
|
4.945235359296727046734888, |
|
133
|
|
|
201.8112620856775083915565, |
|
134
|
|
|
2290.838373831346393026739, |
|
135
|
|
|
11319.67205903380828685045, |
|
136
|
|
|
28557.24635671635335736389, |
|
137
|
|
|
38484.96228443793359990269, |
|
138
|
|
|
26377.48787624195437963534, |
|
139
|
|
|
7225.813979700288197698961, |
|
140
|
|
|
]; |
|
141
|
|
|
|
|
142
|
|
|
private const LG_P2 = [ |
|
143
|
|
|
4.974607845568932035012064, |
|
144
|
|
|
542.4138599891070494101986, |
|
145
|
|
|
15506.93864978364947665077, |
|
146
|
|
|
184793.2904445632425417223, |
|
147
|
|
|
1088204.76946882876749847, |
|
148
|
|
|
3338152.967987029735917223, |
|
149
|
|
|
5106661.678927352456275255, |
|
150
|
|
|
3074109.054850539556250927, |
|
151
|
|
|
]; |
|
152
|
|
|
|
|
153
|
|
|
private const LG_P4 = [ |
|
154
|
|
|
14745.02166059939948905062, |
|
155
|
|
|
2426813.369486704502836312, |
|
156
|
|
|
121475557.4045093227939592, |
|
157
|
|
|
2663432449.630976949898078, |
|
158
|
|
|
29403789566.34553899906876, |
|
159
|
|
|
170266573776.5398868392998, |
|
160
|
|
|
492612579337.743088758812, |
|
161
|
|
|
560625185622.3951465078242, |
|
162
|
|
|
]; |
|
163
|
|
|
|
|
164
|
|
|
private const LG_Q1 = [ |
|
165
|
|
|
67.48212550303777196073036, |
|
166
|
|
|
1113.332393857199323513008, |
|
167
|
|
|
7738.757056935398733233834, |
|
168
|
|
|
27639.87074403340708898585, |
|
169
|
|
|
54993.10206226157329794414, |
|
170
|
|
|
61611.22180066002127833352, |
|
171
|
|
|
36351.27591501940507276287, |
|
172
|
|
|
8785.536302431013170870835, |
|
173
|
|
|
]; |
|
174
|
|
|
|
|
175
|
|
|
private const LG_Q2 = [ |
|
176
|
|
|
183.0328399370592604055942, |
|
177
|
|
|
7765.049321445005871323047, |
|
178
|
|
|
133190.3827966074194402448, |
|
179
|
|
|
1136705.821321969608938755, |
|
180
|
|
|
5267964.117437946917577538, |
|
181
|
|
|
13467014.54311101692290052, |
|
182
|
|
|
17827365.30353274213975932, |
|
183
|
|
|
9533095.591844353613395747, |
|
184
|
|
|
]; |
|
185
|
|
|
|
|
186
|
|
|
private const LG_Q4 = [ |
|
187
|
|
|
2690.530175870899333379843, |
|
188
|
|
|
639388.5654300092398984238, |
|
189
|
|
|
41355999.30241388052042842, |
|
190
|
|
|
1120872109.61614794137657, |
|
191
|
|
|
14886137286.78813811542398, |
|
192
|
|
|
101680358627.2438228077304, |
|
193
|
|
|
341747634550.7377132798597, |
|
194
|
|
|
446315818741.9713286462081, |
|
195
|
|
|
]; |
|
196
|
|
|
|
|
197
|
|
|
private const LG_C = [ |
|
198
|
|
|
-0.001910444077728, |
|
199
|
|
|
8.4171387781295e-4, |
|
200
|
|
|
-5.952379913043012e-4, |
|
201
|
|
|
7.93650793500350248e-4, |
|
202
|
|
|
-0.002777777777777681622553, |
|
203
|
|
|
0.08333333333333333331554247, |
|
204
|
|
|
0.0057083835261, |
|
205
|
|
|
]; |
|
206
|
|
|
|
|
207
|
|
|
// Rough estimate of the fourth root of logGamma_xBig |
|
208
|
|
|
private const LG_FRTBIG = 2.25e76; |
|
209
|
|
|
|
|
210
|
|
|
private const PNT68 = 0.6796875; |
|
211
|
|
|
|
|
212
|
|
|
// Function cache for logGamma |
|
213
|
|
|
|
|
214
|
|
|
private static float $logGammaCacheResult = 0.0; |
|
215
|
|
|
|
|
216
|
|
|
private static float $logGammaCacheX = 0.0; |
|
217
|
|
|
|
|
218
|
|
|
/** |
|
219
|
|
|
* logGamma function. |
|
220
|
|
|
* |
|
221
|
|
|
* Original author was Jaco van Kooten. Ported to PHP by Paul Meagher. |
|
222
|
|
|
* |
|
223
|
|
|
* The natural logarithm of the gamma function. <br /> |
|
224
|
|
|
* Based on public domain NETLIB (Fortran) code by W. J. Cody and L. Stoltz <br /> |
|
225
|
|
|
* Applied Mathematics Division <br /> |
|
226
|
|
|
* Argonne National Laboratory <br /> |
|
227
|
|
|
* Argonne, IL 60439 <br /> |
|
228
|
|
|
* <p> |
|
229
|
|
|
* References: |
|
230
|
|
|
* <ol> |
|
231
|
|
|
* <li>W. J. Cody and K. E. Hillstrom, 'Chebyshev Approximations for the Natural |
|
232
|
|
|
* Logarithm of the Gamma Function,' Math. Comp. 21, 1967, pp. 198-203.</li> |
|
233
|
|
|
* <li>K. E. Hillstrom, ANL/AMD Program ANLC366S, DGAMMA/DLGAMA, May, 1969.</li> |
|
234
|
|
|
* <li>Hart, Et. Al., Computer Approximations, Wiley and sons, New York, 1968.</li> |
|
235
|
|
|
* </ol> |
|
236
|
|
|
* </p> |
|
237
|
|
|
* <p> |
|
238
|
|
|
* From the original documentation: |
|
239
|
|
|
* </p> |
|
240
|
|
|
* <p> |
|
241
|
|
|
* This routine calculates the LOG(GAMMA) function for a positive real argument X. |
|
242
|
|
|
* Computation is based on an algorithm outlined in references 1 and 2. |
|
243
|
|
|
* The program uses rational functions that theoretically approximate LOG(GAMMA) |
|
244
|
|
|
* to at least 18 significant decimal digits. The approximation for X > 12 is from |
|
245
|
|
|
* reference 3, while approximations for X < 12.0 are similar to those in reference |
|
246
|
|
|
* 1, but are unpublished. The accuracy achieved depends on the arithmetic system, |
|
247
|
|
|
* the compiler, the intrinsic functions, and proper selection of the |
|
248
|
|
|
* machine-dependent constants. |
|
249
|
|
|
* </p> |
|
250
|
|
|
* <p> |
|
251
|
|
|
* Error returns: <br /> |
|
252
|
|
|
* The program returns the value XINF for X .LE. 0.0 or when overflow would occur. |
|
253
|
|
|
* The computation is believed to be free of underflow and overflow. |
|
254
|
|
|
* </p> |
|
255
|
|
|
* |
|
256
|
|
|
* @version 1.1 |
|
257
|
|
|
* |
|
258
|
|
|
* @author Jaco van Kooten |
|
259
|
|
|
* |
|
260
|
|
|
* @return float MAX_VALUE for x < 0.0 or when overflow would occur, i.e. x > 2.55E305 |
|
261
|
|
|
*/ |
|
262
|
20 |
|
public static function logGamma(float $x): float |
|
263
|
|
|
{ |
|
264
|
20 |
|
if ($x == self::$logGammaCacheX) { |
|
265
|
1 |
|
return self::$logGammaCacheResult; |
|
266
|
|
|
} |
|
267
|
|
|
|
|
268
|
20 |
|
$y = $x; |
|
269
|
20 |
|
if ($y > 0.0 && $y <= self::LOG_GAMMA_X_MAX_VALUE) { |
|
270
|
20 |
|
if ($y <= self::EPS) { |
|
271
|
|
|
$res = -log($y); |
|
272
|
20 |
|
} elseif ($y <= 1.5) { |
|
273
|
4 |
|
$res = self::logGamma1($y); |
|
274
|
17 |
|
} elseif ($y <= 4.0) { |
|
275
|
8 |
|
$res = self::logGamma2($y); |
|
276
|
16 |
|
} elseif ($y <= 12.0) { |
|
277
|
16 |
|
$res = self::logGamma3($y); |
|
278
|
|
|
} else { |
|
279
|
9 |
|
$res = self::logGamma4($y); |
|
280
|
|
|
} |
|
281
|
|
|
} else { |
|
282
|
|
|
// -------------------------- |
|
283
|
|
|
// Return for bad arguments |
|
284
|
|
|
// -------------------------- |
|
285
|
|
|
$res = self::MAX_VALUE; |
|
286
|
|
|
} |
|
287
|
|
|
|
|
288
|
|
|
// ------------------------------ |
|
289
|
|
|
// Final adjustments and return |
|
290
|
|
|
// ------------------------------ |
|
291
|
20 |
|
self::$logGammaCacheX = $x; |
|
292
|
20 |
|
self::$logGammaCacheResult = $res; |
|
293
|
|
|
|
|
294
|
20 |
|
return $res; |
|
295
|
|
|
} |
|
296
|
|
|
|
|
297
|
4 |
|
private static function logGamma1(float $y): float |
|
298
|
|
|
{ |
|
299
|
|
|
// --------------------- |
|
300
|
|
|
// EPS .LT. X .LE. 1.5 |
|
301
|
|
|
// --------------------- |
|
302
|
4 |
|
if ($y < self::PNT68) { |
|
303
|
3 |
|
$corr = -log($y); |
|
304
|
3 |
|
$xm1 = $y; |
|
305
|
|
|
} else { |
|
306
|
3 |
|
$corr = 0.0; |
|
307
|
3 |
|
$xm1 = $y - 1.0; |
|
308
|
|
|
} |
|
309
|
|
|
|
|
310
|
4 |
|
$xden = 1.0; |
|
311
|
4 |
|
$xnum = 0.0; |
|
312
|
4 |
|
if ($y <= 0.5 || $y >= self::PNT68) { |
|
313
|
4 |
|
for ($i = 0; $i < 8; ++$i) { |
|
314
|
4 |
|
$xnum = $xnum * $xm1 + self::LG_P1[$i]; |
|
315
|
4 |
|
$xden = $xden * $xm1 + self::LG_Q1[$i]; |
|
316
|
|
|
} |
|
317
|
|
|
|
|
318
|
4 |
|
return $corr + $xm1 * (self::LG_D1 + $xm1 * ($xnum / $xden)); |
|
319
|
|
|
} |
|
320
|
|
|
|
|
321
|
|
|
$xm2 = $y - 1.0; |
|
322
|
|
|
for ($i = 0; $i < 8; ++$i) { |
|
323
|
|
|
$xnum = $xnum * $xm2 + self::LG_P2[$i]; |
|
324
|
|
|
$xden = $xden * $xm2 + self::LG_Q2[$i]; |
|
325
|
|
|
} |
|
326
|
|
|
|
|
327
|
|
|
return $corr + $xm2 * (self::LG_D2 + $xm2 * ($xnum / $xden)); |
|
328
|
|
|
} |
|
329
|
|
|
|
|
330
|
8 |
|
private static function logGamma2(float $y): float |
|
331
|
|
|
{ |
|
332
|
|
|
// --------------------- |
|
333
|
|
|
// 1.5 .LT. X .LE. 4.0 |
|
334
|
|
|
// --------------------- |
|
335
|
8 |
|
$xm2 = $y - 2.0; |
|
336
|
8 |
|
$xden = 1.0; |
|
337
|
8 |
|
$xnum = 0.0; |
|
338
|
8 |
|
for ($i = 0; $i < 8; ++$i) { |
|
339
|
8 |
|
$xnum = $xnum * $xm2 + self::LG_P2[$i]; |
|
340
|
8 |
|
$xden = $xden * $xm2 + self::LG_Q2[$i]; |
|
341
|
|
|
} |
|
342
|
|
|
|
|
343
|
8 |
|
return $xm2 * (self::LG_D2 + $xm2 * ($xnum / $xden)); |
|
344
|
|
|
} |
|
345
|
|
|
|
|
346
|
16 |
|
protected static function logGamma3(float $y): float |
|
347
|
|
|
{ |
|
348
|
|
|
// ---------------------- |
|
349
|
|
|
// 4.0 .LT. X .LE. 12.0 |
|
350
|
|
|
// ---------------------- |
|
351
|
16 |
|
$xm4 = $y - 4.0; |
|
352
|
16 |
|
$xden = -1.0; |
|
353
|
16 |
|
$xnum = 0.0; |
|
354
|
16 |
|
for ($i = 0; $i < 8; ++$i) { |
|
355
|
16 |
|
$xnum = $xnum * $xm4 + self::LG_P4[$i]; |
|
356
|
16 |
|
$xden = $xden * $xm4 + self::LG_Q4[$i]; |
|
357
|
|
|
} |
|
358
|
|
|
|
|
359
|
16 |
|
return self::LG_D4 + $xm4 * ($xnum / $xden); |
|
360
|
|
|
} |
|
361
|
|
|
|
|
362
|
9 |
|
protected static function logGamma4(float $y): float |
|
363
|
|
|
{ |
|
364
|
|
|
// --------------------------------- |
|
365
|
|
|
// Evaluate for argument .GE. 12.0 |
|
366
|
|
|
// --------------------------------- |
|
367
|
9 |
|
$res = 0.0; |
|
368
|
9 |
|
if ($y <= self::LG_FRTBIG) { |
|
369
|
9 |
|
$res = self::LG_C[6]; |
|
370
|
9 |
|
$ysq = $y * $y; |
|
371
|
9 |
|
for ($i = 0; $i < 6; ++$i) { |
|
372
|
9 |
|
$res = $res / $ysq + self::LG_C[$i]; |
|
373
|
|
|
} |
|
374
|
9 |
|
$res /= $y; |
|
375
|
9 |
|
$corr = log($y); |
|
376
|
9 |
|
$res = $res + log(self::SQRT2PI) - 0.5 * $corr; |
|
377
|
9 |
|
$res += $y * ($corr - 1.0); |
|
378
|
|
|
} |
|
379
|
|
|
|
|
380
|
9 |
|
return $res; |
|
381
|
|
|
} |
|
382
|
|
|
} |
|
383
|
|
|
|
PHPOffice /
PhpSpreadsheet
