ILIAS  release_5-4 Revision v5.4.26-12-gabc799a52e6
GammaBase.php
Go to the documentation of this file.
1<?php
2
4
6
7abstract class GammaBase
8{
9 private const LOG_GAMMA_X_MAX_VALUE = 2.55e305;
10
11 private const EPS = 2.22e-16;
12
13 private const MAX_VALUE = 1.2e308;
14
15 private const SQRT2PI = 2.5066282746310005024157652848110452530069867406099;
16
17 private const MAX_ITERATIONS = 256;
18
19 protected static function calculateDistribution(float $value, float $a, float $b, bool $cumulative)
20 {
21 if ($cumulative) {
22 return self::incompleteGamma($a, $value / $b) / self::gammaValue($a);
23 }
24
25 return (1 / ($b ** $a * self::gammaValue($a))) * $value ** ($a - 1) * exp(0 - ($value / $b));
26 }
27
28 protected static function calculateInverse(float $probability, float $alpha, float $beta)
29 {
30 $xLo = 0;
31 $xHi = $alpha * $beta * 5;
32
33 $dx = 1024;
34 $x = $xNew = 1;
35 $i = 0;
36
37 while ((abs($dx) > Functions::PRECISION) && (++$i <= self::MAX_ITERATIONS)) {
38 // Apply Newton-Raphson step
39 $result = self::calculateDistribution($x, $alpha, $beta, true);
40 $error = $result - $probability;
41
42 if ($error == 0.0) {
43 $dx = 0;
44 } elseif ($error < 0.0) {
45 $xLo = $x;
46 } else {
47 $xHi = $x;
48 }
49
50 $pdf = self::calculateDistribution($x, $alpha, $beta, false);
51 // Avoid division by zero
52 if ($pdf !== 0.0) {
53 $dx = $error / $pdf;
54 $xNew = $x - $dx;
55 }
56
57 // If the NR fails to converge (which for example may be the
58 // case if the initial guess is too rough) we apply a bisection
59 // step to determine a more narrow interval around the root.
60 if (($xNew < $xLo) || ($xNew > $xHi) || ($pdf == 0.0)) {
61 $xNew = ($xLo + $xHi) / 2;
62 $dx = $xNew - $x;
63 }
64 $x = $xNew;
65 }
66
67 if ($i === self::MAX_ITERATIONS) {
68 return Functions::NA();
69 }
70
71 return $x;
72 }
73
74 //
75 // Implementation of the incomplete Gamma function
76 //
77 public static function incompleteGamma(float $a, float $x): float
78 {
79 static $max = 32;
80 $summer = 0;
81 for ($n = 0; $n <= $max; ++$n) {
82 $divisor = $a;
83 for ($i = 1; $i <= $n; ++$i) {
84 $divisor *= ($a + $i);
85 }
86 $summer += ($x ** $n / $divisor);
87 }
88
89 return $x ** $a * exp(0 - $x) * $summer;
90 }
91
92 //
93 // Implementation of the Gamma function
94 //
95 public static function gammaValue(float $value): float
96 {
97 if ($value == 0.0) {
98 return 0;
99 }
100
101 static $p0 = 1.000000000190015;
102 static $p = [
103 1 => 76.18009172947146,
104 2 => -86.50532032941677,
105 3 => 24.01409824083091,
106 4 => -1.231739572450155,
107 5 => 1.208650973866179e-3,
108 6 => -5.395239384953e-6,
109 ];
110
111 $y = $x = $value;
112 $tmp = $x + 5.5;
113 $tmp -= ($x + 0.5) * log($tmp);
114
115 $summer = $p0;
116 for ($j = 1; $j <= 6; ++$j) {
117 $summer += ($p[$j] / ++$y);
118 }
119
120 return exp(0 - $tmp + log(self::SQRT2PI * $summer / $x));
121 }
122
168 // Log Gamma related constants
169 private const LG_D1 = -0.5772156649015328605195174;
170
171 private const LG_D2 = 0.4227843350984671393993777;
172
173 private const LG_D4 = 1.791759469228055000094023;
174
175 private const LG_P1 = [
176 4.945235359296727046734888,
177 201.8112620856775083915565,
178 2290.838373831346393026739,
179 11319.67205903380828685045,
180 28557.24635671635335736389,
181 38484.96228443793359990269,
182 26377.48787624195437963534,
183 7225.813979700288197698961,
184 ];
185
186 private const LG_P2 = [
187 4.974607845568932035012064,
188 542.4138599891070494101986,
189 15506.93864978364947665077,
190 184793.2904445632425417223,
191 1088204.76946882876749847,
192 3338152.967987029735917223,
193 5106661.678927352456275255,
194 3074109.054850539556250927,
195 ];
196
197 private const LG_P4 = [
198 14745.02166059939948905062,
199 2426813.369486704502836312,
200 121475557.4045093227939592,
201 2663432449.630976949898078,
202 29403789566.34553899906876,
203 170266573776.5398868392998,
204 492612579337.743088758812,
205 560625185622.3951465078242,
206 ];
207
208 private const LG_Q1 = [
209 67.48212550303777196073036,
210 1113.332393857199323513008,
211 7738.757056935398733233834,
212 27639.87074403340708898585,
213 54993.10206226157329794414,
214 61611.22180066002127833352,
215 36351.27591501940507276287,
216 8785.536302431013170870835,
217 ];
218
219 private const LG_Q2 = [
220 183.0328399370592604055942,
221 7765.049321445005871323047,
222 133190.3827966074194402448,
223 1136705.821321969608938755,
224 5267964.117437946917577538,
225 13467014.54311101692290052,
226 17827365.30353274213975932,
227 9533095.591844353613395747,
228 ];
229
230 private const LG_Q4 = [
231 2690.530175870899333379843,
232 639388.5654300092398984238,
233 41355999.30241388052042842,
234 1120872109.61614794137657,
235 14886137286.78813811542398,
236 101680358627.2438228077304,
237 341747634550.7377132798597,
238 446315818741.9713286462081,
239 ];
240
241 private const LG_C = [
242 -0.001910444077728,
243 8.4171387781295e-4,
244 -5.952379913043012e-4,
245 7.93650793500350248e-4,
246 -0.002777777777777681622553,
247 0.08333333333333333331554247,
248 0.0057083835261,
249 ];
250
251 // Rough estimate of the fourth root of logGamma_xBig
252 private const LG_FRTBIG = 2.25e76;
253
254 private const PNT68 = 0.6796875;
255
256 // Function cache for logGamma
257 private static $logGammaCacheResult = 0.0;
258
259 private static $logGammaCacheX = 0.0;
260
261 public static function logGamma(float $x): float
262 {
263 if ($x == self::$logGammaCacheX) {
265 }
266
267 $y = $x;
268 if ($y > 0.0 && $y <= self::LOG_GAMMA_X_MAX_VALUE) {
269 if ($y <= self::EPS) {
270 $res = -log($y);
271 } elseif ($y <= 1.5) {
273 } elseif ($y <= 4.0) {
275 } elseif ($y <= 12.0) {
277 } else {
279 }
280 } else {
281 // --------------------------
282 // Return for bad arguments
283 // --------------------------
285 }
286
287 // ------------------------------
288 // Final adjustments and return
289 // ------------------------------
290 self::$logGammaCacheX = $x;
291 self::$logGammaCacheResult = $res;
292
293 return $res;
294 }
295
296 private static function logGamma1(float $y)
297 {
298 // ---------------------
299 // EPS .LT. X .LE. 1.5
300 // ---------------------
301 if ($y < self::PNT68) {
302 $corr = -log($y);
303 $xm1 = $y;
304 } else {
305 $corr = 0.0;
306 $xm1 = $y - 1.0;
307 }
308
309 $xden = 1.0;
310 $xnum = 0.0;
311 if ($y <= 0.5 || $y >= self::PNT68) {
312 for ($i = 0; $i < 8; ++$i) {
313 $xnum = $xnum * $xm1 + self::LG_P1[$i];
314 $xden = $xden * $xm1 + self::LG_Q1[$i];
315 }
316
317 return $corr + $xm1 * (self::LG_D1 + $xm1 * ($xnum / $xden));
318 }
319
320 $xm2 = $y - 1.0;
321 for ($i = 0; $i < 8; ++$i) {
322 $xnum = $xnum * $xm2 + self::LG_P2[$i];
323 $xden = $xden * $xm2 + self::LG_Q2[$i];
324 }
325
326 return $corr + $xm2 * (self::LG_D2 + $xm2 * ($xnum / $xden));
327 }
328
329 private static function logGamma2(float $y)
330 {
331 // ---------------------
332 // 1.5 .LT. X .LE. 4.0
333 // ---------------------
334 $xm2 = $y - 2.0;
335 $xden = 1.0;
336 $xnum = 0.0;
337 for ($i = 0; $i < 8; ++$i) {
338 $xnum = $xnum * $xm2 + self::LG_P2[$i];
339 $xden = $xden * $xm2 + self::LG_Q2[$i];
340 }
341
342 return $xm2 * (self::LG_D2 + $xm2 * ($xnum / $xden));
343 }
344
345 protected static function logGamma3(float $y)
346 {
347 // ----------------------
348 // 4.0 .LT. X .LE. 12.0
349 // ----------------------
350 $xm4 = $y - 4.0;
351 $xden = -1.0;
352 $xnum = 0.0;
353 for ($i = 0; $i < 8; ++$i) {
354 $xnum = $xnum * $xm4 + self::LG_P4[$i];
355 $xden = $xden * $xm4 + self::LG_Q4[$i];
356 }
357
358 return self::LG_D4 + $xm4 * ($xnum / $xden);
359 }
360
361 protected static function logGamma4(float $y)
362 {
363 // ---------------------------------
364 // Evaluate for argument .GE. 12.0
365 // ---------------------------------
366 $res = 0.0;
367 if ($y <= self::LG_FRTBIG) {
368 $res = self::LG_C[6];
369 $ysq = $y * $y;
370 for ($i = 0; $i < 6; ++$i) {
371 $res = $res / $ysq + self::LG_C[$i];
372 }
373 $res /= $y;
374 $corr = log($y);
375 $res = $res + log(self::SQRT2PI) - 0.5 * $corr;
376 $res += $y * ($corr - 1.0);
377 }
378
379 return $res;
380 }
381}
$result
$n
Definition: RandomTest.php:85
An exception for terminatinating execution or to throw for unit testing.
static calculateDistribution(float $value, float $a, float $b, bool $cumulative)
Definition: GammaBase.php:19
static calculateInverse(float $probability, float $alpha, float $beta)
Definition: GammaBase.php:28
$x
Definition: complexTest.php:9
$i
Definition: disco.tpl.php:19
$pdf
Definition: example_001.php:31
$y
Definition: example_007.php:83
foreach($_POST as $key=> $value) $res