[Summary view]

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

   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 ChiSquared
  11  {
  12      use ArrayEnabled;
  13  
  14      private const MAX_ITERATIONS = 256;
  15  
  16      private const EPS = 2.22e-16;
  17  
  18      /**
  19       * CHIDIST.
  20       *
  21       * Returns the one-tailed probability of the chi-squared distribution.
  22       *
  23       * @param mixed $value Float value for which we want the probability
  24       *                      Or can be an array of values
  25       * @param mixed $degrees Integer degrees of freedom
  26       *                      Or can be an array of values
  27       *
  28       * @return array|float|string
  29       *         If an array of numbers is passed as an argument, then the returned result will also be an array
  30       *            with the same dimensions
  31       */
  32      public static function distributionRightTail($value, $degrees)
  33      {
  34          if (is_array($value) || is_array($degrees)) {
  35              return self::evaluateArrayArguments([self::class, __FUNCTION__], $value, $degrees);
  36          }
  37  
  38          try {
  39              $value = DistributionValidations::validateFloat($value);
  40              $degrees = DistributionValidations::validateInt($degrees);
  41          } catch (Exception $e) {
  42              return $e->getMessage();
  43          }
  44  
  45          if ($degrees < 1) {
  46              return ExcelError::NAN();
  47          }
  48          if ($value < 0) {
  49              if (Functions::getCompatibilityMode() == Functions::COMPATIBILITY_GNUMERIC) {
  50                  return 1;
  51              }
  52  
  53              return ExcelError::NAN();
  54          }
  55  
  56          return 1 - (Gamma::incompleteGamma($degrees / 2, $value / 2) / Gamma::gammaValue($degrees / 2));
  57      }
  58  
  59      /**
  60       * CHIDIST.
  61       *
  62       * Returns the one-tailed probability of the chi-squared distribution.
  63       *
  64       * @param mixed $value Float value for which we want the probability
  65       *                      Or can be an array of values
  66       * @param mixed $degrees Integer degrees of freedom
  67       *                      Or can be an array of values
  68       * @param mixed $cumulative Boolean value indicating if we want the cdf (true) or the pdf (false)
  69       *                      Or can be an array of values
  70       *
  71       * @return array|float|string
  72       *         If an array of numbers is passed as an argument, then the returned result will also be an array
  73       *            with the same dimensions
  74       */
  75      public static function distributionLeftTail($value, $degrees, $cumulative)
  76      {
  77          if (is_array($value) || is_array($degrees) || is_array($cumulative)) {
  78              return self::evaluateArrayArguments([self::class, __FUNCTION__], $value, $degrees, $cumulative);
  79          }
  80  
  81          try {
  82              $value = DistributionValidations::validateFloat($value);
  83              $degrees = DistributionValidations::validateInt($degrees);
  84              $cumulative = DistributionValidations::validateBool($cumulative);
  85          } catch (Exception $e) {
  86              return $e->getMessage();
  87          }
  88  
  89          if ($degrees < 1) {
  90              return ExcelError::NAN();
  91          }
  92          if ($value < 0) {
  93              if (Functions::getCompatibilityMode() == Functions::COMPATIBILITY_GNUMERIC) {
  94                  return 1;
  95              }
  96  
  97              return ExcelError::NAN();
  98          }
  99  
 100          if ($cumulative === true) {
 101              return 1 - self::distributionRightTail($value, $degrees);
 102          }
 103  
 104          return ($value ** (($degrees / 2) - 1) * exp(-$value / 2)) /
 105              ((2 ** ($degrees / 2)) * Gamma::gammaValue($degrees / 2));
 106      }
 107  
 108      /**
 109       * CHIINV.
 110       *
 111       * Returns the inverse of the right-tailed probability of the chi-squared distribution.
 112       *
 113       * @param mixed $probability Float probability at which you want to evaluate the distribution
 114       *                      Or can be an array of values
 115       * @param mixed $degrees Integer degrees of freedom
 116       *                      Or can be an array of values
 117       *
 118       * @return array|float|string
 119       *         If an array of numbers is passed as an argument, then the returned result will also be an array
 120       *            with the same dimensions
 121       */
 122      public static function inverseRightTail($probability, $degrees)
 123      {
 124          if (is_array($probability) || is_array($degrees)) {
 125              return self::evaluateArrayArguments([self::class, __FUNCTION__], $probability, $degrees);
 126          }
 127  
 128          try {
 129              $probability = DistributionValidations::validateProbability($probability);
 130              $degrees = DistributionValidations::validateInt($degrees);
 131          } catch (Exception $e) {
 132              return $e->getMessage();
 133          }
 134  
 135          if ($degrees < 1) {
 136              return ExcelError::NAN();
 137          }
 138  
 139          $callback = function ($value) use ($degrees) {
 140              return 1 - (Gamma::incompleteGamma($degrees / 2, $value / 2)
 141                      / Gamma::gammaValue($degrees / 2));
 142          };
 143  
 144          $newtonRaphson = new NewtonRaphson($callback);
 145  
 146          return $newtonRaphson->execute($probability);
 147      }
 148  
 149      /**
 150       * CHIINV.
 151       *
 152       * Returns the inverse of the left-tailed probability of the chi-squared distribution.
 153       *
 154       * @param mixed $probability Float probability at which you want to evaluate the distribution
 155       *                      Or can be an array of values
 156       * @param mixed $degrees Integer degrees of freedom
 157       *                      Or can be an array of values
 158       *
 159       * @return array|float|string
 160       *         If an array of numbers is passed as an argument, then the returned result will also be an array
 161       *            with the same dimensions
 162       */
 163      public static function inverseLeftTail($probability, $degrees)
 164      {
 165          if (is_array($probability) || is_array($degrees)) {
 166              return self::evaluateArrayArguments([self::class, __FUNCTION__], $probability, $degrees);
 167          }
 168  
 169          try {
 170              $probability = DistributionValidations::validateProbability($probability);
 171              $degrees = DistributionValidations::validateInt($degrees);
 172          } catch (Exception $e) {
 173              return $e->getMessage();
 174          }
 175  
 176          if ($degrees < 1) {
 177              return ExcelError::NAN();
 178          }
 179  
 180          return self::inverseLeftTailCalculation($probability, $degrees);
 181      }
 182  
 183      /**
 184       * CHITEST.
 185       *
 186       * Uses the chi-square test to calculate the probability that the differences between two supplied data sets
 187       *      (of observed and expected frequencies), are likely to be simply due to sampling error,
 188       *      or if they are likely to be real.
 189       *
 190       * @param mixed $actual an array of observed frequencies
 191       * @param mixed $expected an array of expected frequencies
 192       *
 193       * @return float|string
 194       */
 195      public static function test($actual, $expected)
 196      {
 197          $rows = count($actual);
 198          $actual = Functions::flattenArray($actual);
 199          $expected = Functions::flattenArray($expected);
 200          $columns = count($actual) / $rows;
 201  
 202          $countActuals = count($actual);
 203          $countExpected = count($expected);
 204          if ($countActuals !== $countExpected || $countActuals === 1) {
 205              return ExcelError::NAN();
 206          }
 207  
 208          $result = 0.0;
 209          for ($i = 0; $i < $countActuals; ++$i) {
 210              if ($expected[$i] == 0.0) {
 211                  return ExcelError::DIV0();
 212              } elseif ($expected[$i] < 0.0) {
 213                  return ExcelError::NAN();
 214              }
 215              $result += (($actual[$i] - $expected[$i]) ** 2) / $expected[$i];
 216          }
 217  
 218          $degrees = self::degrees($rows, $columns);
 219  
 220          $result = Functions::scalar(self::distributionRightTail($result, $degrees));
 221  
 222          return $result;
 223      }
 224  
 225      protected static function degrees(int $rows, int $columns): int
 226      {
 227          if ($rows === 1) {
 228              return $columns - 1;
 229          } elseif ($columns === 1) {
 230              return $rows - 1;
 231          }
 232  
 233          return ($columns - 1) * ($rows - 1);
 234      }
 235  
 236      private static function inverseLeftTailCalculation(float $probability, int $degrees): float
 237      {
 238          // bracket the root
 239          $min = 0;
 240          $sd = sqrt(2.0 * $degrees);
 241          $max = 2 * $sd;
 242          $s = -1;
 243  
 244          while ($s * self::pchisq($max, $degrees) > $probability * $s) {
 245              $min = $max;
 246              $max += 2 * $sd;
 247          }
 248  
 249          // Find root using bisection
 250          $chi2 = 0.5 * ($min + $max);
 251  
 252          while (($max - $min) > self::EPS * $chi2) {
 253              if ($s * self::pchisq($chi2, $degrees) > $probability * $s) {
 254                  $min = $chi2;
 255              } else {
 256                  $max = $chi2;
 257              }
 258              $chi2 = 0.5 * ($min + $max);
 259          }
 260  
 261          return $chi2;
 262      }
 263  
 264      private static function pchisq($chi2, $degrees)
 265      {
 266          return self::gammp($degrees, 0.5 * $chi2);
 267      }
 268  
 269      private static function gammp($n, $x)
 270      {
 271          if ($x < 0.5 * $n + 1) {
 272              return self::gser($n, $x);
 273          }
 274  
 275          return 1 - self::gcf($n, $x);
 276      }
 277  
 278      // Return the incomplete gamma function P(n/2,x) evaluated by
 279      // series representation. Algorithm from numerical recipe.
 280      // Assume that n is a positive integer and x>0, won't check arguments.
 281      // Relative error controlled by the eps parameter
 282      private static function gser($n, $x)
 283      {
 284          /** @var float */
 285          $gln = Gamma::ln($n / 2);
 286          $a = 0.5 * $n;
 287          $ap = $a;
 288          $sum = 1.0 / $a;
 289          $del = $sum;
 290          for ($i = 1; $i < 101; ++$i) {
 291              ++$ap;
 292              $del = $del * $x / $ap;
 293              $sum += $del;
 294              if ($del < $sum * self::EPS) {
 295                  break;
 296              }
 297          }
 298  
 299          return $sum * exp(-$x + $a * log($x) - $gln);
 300      }
 301  
 302      // Return the incomplete gamma function Q(n/2,x) evaluated by
 303      // its continued fraction representation. Algorithm from numerical recipe.
 304      // Assume that n is a postive integer and x>0, won't check arguments.
 305      // Relative error controlled by the eps parameter
 306      private static function gcf($n, $x)
 307      {
 308          /** @var float */
 309          $gln = Gamma::ln($n / 2);
 310          $a = 0.5 * $n;
 311          $b = $x + 1 - $a;
 312          $fpmin = 1.e-300;
 313          $c = 1 / $fpmin;
 314          $d = 1 / $b;
 315          $h = $d;
 316          for ($i = 1; $i < 101; ++$i) {
 317              $an = -$i * ($i - $a);
 318              $b += 2;
 319              $d = $an * $d + $b;
 320              if (abs($d) < $fpmin) {
 321                  $d = $fpmin;
 322              }
 323              $c = $b + $an / $c;
 324              if (abs($c) < $fpmin) {
 325                  $c = $fpmin;
 326              }
 327              $d = 1 / $d;
 328              $del = $d * $c;
 329              $h = $h * $del;
 330              if (abs($del - 1) < self::EPS) {
 331                  break;
 332              }
 333          }
 334  
 335          return $h * exp(-$x + $a * log($x) - $gln);
 336      }
 337  }