[Summary view]

Differences Between: [Versions 400 and 403] [Versions 401 and 403]

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