[Summary view]

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

   1  <?php
   2  
   3  namespace PhpOffice\PhpSpreadsheet\Calculation\Statistical\Distributions;
   4  
   5  use PhpOffice\PhpSpreadsheet\Calculation\ArrayEnabled;
   6  use PhpOffice\PhpSpreadsheet\Calculation\Exception;
   7  use PhpOffice\PhpSpreadsheet\Calculation\Functions;
   8  use PhpOffice\PhpSpreadsheet\Calculation\Information\ExcelError;
   9  
  10  class Beta
  11  {
  12      use ArrayEnabled;
  13  
  14      private const MAX_ITERATIONS = 256;
  15  
  16      private const LOG_GAMMA_X_MAX_VALUE = 2.55e305;
  17  
  18      private const XMININ = 2.23e-308;
  19  
  20      /**
  21       * BETADIST.
  22       *
  23       * Returns the beta distribution.
  24       *
  25       * @param mixed $value Float value at which you want to evaluate the distribution
  26       *                      Or can be an array of values
  27       * @param mixed $alpha Parameter to the distribution as a float
  28       *                      Or can be an array of values
  29       * @param mixed $beta Parameter to the distribution as a float
  30       *                      Or can be an array of values
  31       * @param mixed $rMin as an float
  32       *                      Or can be an array of values
  33       * @param mixed $rMax as an float
  34       *                      Or can be an array of values
  35       *
  36       * @return array|float|string
  37       *         If an array of numbers is passed as an argument, then the returned result will also be an array
  38       *            with the same dimensions
  39       */
  40      public static function distribution($value, $alpha, $beta, $rMin = 0.0, $rMax = 1.0)
  41      {
  42          if (is_array($value) || is_array($alpha) || is_array($beta) || is_array($rMin) || is_array($rMax)) {
  43              return self::evaluateArrayArguments([self::class, __FUNCTION__], $value, $alpha, $beta, $rMin, $rMax);
  44          }
  45  
  46          $rMin = $rMin ?? 0.0;
  47          $rMax = $rMax ?? 1.0;
  48  
  49          try {
  50              $value = DistributionValidations::validateFloat($value);
  51              $alpha = DistributionValidations::validateFloat($alpha);
  52              $beta = DistributionValidations::validateFloat($beta);
  53              $rMax = DistributionValidations::validateFloat($rMax);
  54              $rMin = DistributionValidations::validateFloat($rMin);
  55          } catch (Exception $e) {
  56              return $e->getMessage();
  57          }
  58  
  59          if ($rMin > $rMax) {
  60              $tmp = $rMin;
  61              $rMin = $rMax;
  62              $rMax = $tmp;
  63          }
  64          if (($value < $rMin) || ($value > $rMax) || ($alpha <= 0) || ($beta <= 0) || ($rMin == $rMax)) {
  65              return ExcelError::NAN();
  66          }
  67  
  68          $value -= $rMin;
  69          $value /= ($rMax - $rMin);
  70  
  71          return self::incompleteBeta($value, $alpha, $beta);
  72      }
  73  
  74      /**
  75       * BETAINV.
  76       *
  77       * Returns the inverse of the Beta distribution.
  78       *
  79       * @param mixed $probability Float probability at which you want to evaluate the distribution
  80       *                      Or can be an array of values
  81       * @param mixed $alpha Parameter to the distribution as a float
  82       *                      Or can be an array of values
  83       * @param mixed $beta Parameter to the distribution as a float
  84       *                      Or can be an array of values
  85       * @param mixed $rMin Minimum value as a float
  86       *                      Or can be an array of values
  87       * @param mixed $rMax Maximum value as a float
  88       *                      Or can be an array of values
  89       *
  90       * @return array|float|string
  91       *         If an array of numbers is passed as an argument, then the returned result will also be an array
  92       *            with the same dimensions
  93       */
  94      public static function inverse($probability, $alpha, $beta, $rMin = 0.0, $rMax = 1.0)
  95      {
  96          if (is_array($probability) || is_array($alpha) || is_array($beta) || is_array($rMin) || is_array($rMax)) {
  97              return self::evaluateArrayArguments([self::class, __FUNCTION__], $probability, $alpha, $beta, $rMin, $rMax);
  98          }
  99  
 100          $rMin = $rMin ?? 0.0;
 101          $rMax = $rMax ?? 1.0;
 102  
 103          try {
 104              $probability = DistributionValidations::validateProbability($probability);
 105              $alpha = DistributionValidations::validateFloat($alpha);
 106              $beta = DistributionValidations::validateFloat($beta);
 107              $rMax = DistributionValidations::validateFloat($rMax);
 108              $rMin = DistributionValidations::validateFloat($rMin);
 109          } catch (Exception $e) {
 110              return $e->getMessage();
 111          }
 112  
 113          if ($rMin > $rMax) {
 114              $tmp = $rMin;
 115              $rMin = $rMax;
 116              $rMax = $tmp;
 117          }
 118          if (($alpha <= 0) || ($beta <= 0) || ($rMin == $rMax) || ($probability <= 0.0)) {
 119              return ExcelError::NAN();
 120          }
 121  
 122          return self::calculateInverse($probability, $alpha, $beta, $rMin, $rMax);
 123      }
 124  
 125      /**
 126       * @return float|string
 127       */
 128      private static function calculateInverse(float $probability, float $alpha, float $beta, float $rMin, float $rMax)
 129      {
 130          $a = 0;
 131          $b = 2;
 132          $guess = ($a + $b) / 2;
 133  
 134          $i = 0;
 135          while ((($b - $a) > Functions::PRECISION) && (++$i <= self::MAX_ITERATIONS)) {
 136              $guess = ($a + $b) / 2;
 137              $result = self::distribution($guess, $alpha, $beta);
 138              if (($result === $probability) || ($result === 0.0)) {
 139                  $b = $a;
 140              } elseif ($result > $probability) {
 141                  $b = $guess;
 142              } else {
 143                  $a = $guess;
 144              }
 145          }
 146  
 147          if ($i === self::MAX_ITERATIONS) {
 148              return ExcelError::NA();
 149          }
 150  
 151          return round($rMin + $guess * ($rMax - $rMin), 12);
 152      }
 153  
 154      /**
 155       * Incomplete beta function.
 156       *
 157       * @author Jaco van Kooten
 158       * @author Paul Meagher
 159       *
 160       * The computation is based on formulas from Numerical Recipes, Chapter 6.4 (W.H. Press et al, 1992).
 161       *
 162       * @param float $x require 0<=x<=1
 163       * @param float $p require p>0
 164       * @param float $q require q>0
 165       *
 166       * @return float 0 if x<0, p<=0, q<=0 or p+q>2.55E305 and 1 if x>1 to avoid errors and over/underflow
 167       */
 168      public static function incompleteBeta(float $x, float $p, float $q): float
 169      {
 170          if ($x <= 0.0) {
 171              return 0.0;
 172          } elseif ($x >= 1.0) {
 173              return 1.0;
 174          } elseif (($p <= 0.0) || ($q <= 0.0) || (($p + $q) > self::LOG_GAMMA_X_MAX_VALUE)) {
 175              return 0.0;
 176          }
 177  
 178          $beta_gam = exp((0 - self::logBeta($p, $q)) + $p * log($x) + $q * log(1.0 - $x));
 179          if ($x < ($p + 1.0) / ($p + $q + 2.0)) {
 180              return $beta_gam * self::betaFraction($x, $p, $q) / $p;
 181          }
 182  
 183          return 1.0 - ($beta_gam * self::betaFraction(1 - $x, $q, $p) / $q);
 184      }
 185  
 186      // Function cache for logBeta function
 187      /** @var float */
 188      private static $logBetaCacheP = 0.0;
 189  
 190      /** @var float */
 191      private static $logBetaCacheQ = 0.0;
 192  
 193      /** @var float */
 194      private static $logBetaCacheResult = 0.0;
 195  
 196      /**
 197       * The natural logarithm of the beta function.
 198       *
 199       * @param float $p require p>0
 200       * @param float $q require q>0
 201       *
 202       * @return float 0 if p<=0, q<=0 or p+q>2.55E305 to avoid errors and over/underflow
 203       *
 204       * @author Jaco van Kooten
 205       */
 206      private static function logBeta(float $p, float $q): float
 207      {
 208          if ($p != self::$logBetaCacheP || $q != self::$logBetaCacheQ) {
 209              self::$logBetaCacheP = $p;
 210              self::$logBetaCacheQ = $q;
 211              if (($p <= 0.0) || ($q <= 0.0) || (($p + $q) > self::LOG_GAMMA_X_MAX_VALUE)) {
 212                  self::$logBetaCacheResult = 0.0;
 213              } else {
 214                  self::$logBetaCacheResult = Gamma::logGamma($p) + Gamma::logGamma($q) - Gamma::logGamma($p + $q);
 215              }
 216          }
 217  
 218          return self::$logBetaCacheResult;
 219      }
 220  
 221      /**
 222       * Evaluates of continued fraction part of incomplete beta function.
 223       * Based on an idea from Numerical Recipes (W.H. Press et al, 1992).
 224       *
 225       * @author Jaco van Kooten
 226       */
 227      private static function betaFraction(float $x, float $p, float $q): float
 228      {
 229          $c = 1.0;
 230          $sum_pq = $p + $q;
 231          $p_plus = $p + 1.0;
 232          $p_minus = $p - 1.0;
 233          $h = 1.0 - $sum_pq * $x / $p_plus;
 234          if (abs($h) < self::XMININ) {
 235              $h = self::XMININ;
 236          }
 237          $h = 1.0 / $h;
 238          $frac = $h;
 239          $m = 1;
 240          $delta = 0.0;
 241          while ($m <= self::MAX_ITERATIONS && abs($delta - 1.0) > Functions::PRECISION) {
 242              $m2 = 2 * $m;
 243              // even index for d
 244              $d = $m * ($q - $m) * $x / (($p_minus + $m2) * ($p + $m2));
 245              $h = 1.0 + $d * $h;
 246              if (abs($h) < self::XMININ) {
 247                  $h = self::XMININ;
 248              }
 249              $h = 1.0 / $h;
 250              $c = 1.0 + $d / $c;
 251              if (abs($c) < self::XMININ) {
 252                  $c = self::XMININ;
 253              }
 254              $frac *= $h * $c;
 255              // odd index for d
 256              $d = -($p + $m) * ($sum_pq + $m) * $x / (($p + $m2) * ($p_plus + $m2));
 257              $h = 1.0 + $d * $h;
 258              if (abs($h) < self::XMININ) {
 259                  $h = self::XMININ;
 260              }
 261              $h = 1.0 / $h;
 262              $c = 1.0 + $d / $c;
 263              if (abs($c) < self::XMININ) {
 264                  $c = self::XMININ;
 265              }
 266              $delta = $h * $c;
 267              $frac *= $delta;
 268              ++$m;
 269          }
 270  
 271          return $frac;
 272      }
 273  
 274      /*
 275      private static function betaValue(float $a, float $b): float
 276      {
 277          return (Gamma::gammaValue($a) * Gamma::gammaValue($b)) /
 278              Gamma::gammaValue($a + $b);
 279      }
 280  
 281      private static function regularizedIncompleteBeta(float $value, float $a, float $b): float
 282      {
 283          return self::incompleteBeta($value, $a, $b) / self::betaValue($a, $b);
 284      }
 285      */
 286  }