Skip to content

Commit

Permalink
Chi squared inverse left tailed (#1964)
Browse files Browse the repository at this point in the history
* Implementation of the CHISQ.INV() method for ChiSquared distribution left-tail
  • Loading branch information
Mark Baker authored Mar 28, 2021
1 parent e2ff14f commit e68978f
Show file tree
Hide file tree
Showing 5 changed files with 234 additions and 2 deletions.
2 changes: 1 addition & 1 deletion CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org).

### Added

- Implemented the CHITEST() and CHISQ.DIST() Statistical function.
- Implemented the CHITEST(), CHISQ.DIST() and CHISQ.INV() and equivalent Statistical functions, for both left- and right-tailed distributions.
- Support for ActiveSheet and SelectedCells in the ODS Reader and Writer. [PR #1908](https://github.com/PHPOffice/PhpSpreadsheet/pull/1908)

### Changed
Expand Down
2 changes: 1 addition & 1 deletion src/PhpSpreadsheet/Calculation/Calculation.php
Original file line number Diff line number Diff line change
Expand Up @@ -508,7 +508,7 @@ class Calculation
],
'CHISQ.INV' => [
'category' => Category::CATEGORY_STATISTICAL,
'functionCall' => [Functions::class, 'DUMMY'],
'functionCall' => [Statistical\Distributions\ChiSquared::class, 'inverseLeftTail'],
'argumentCount' => '2',
],
'CHISQ.INV.RT' => [
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,8 @@ class ChiSquared

private const MAX_ITERATIONS = 256;

private const EPS = 2.22e-16;

/**
* CHIDIST.
*
Expand Down Expand Up @@ -127,6 +129,35 @@ public static function inverseRightTail($probability, $degrees)
return $newtonRaphson->execute($probability);
}

/**
* CHIINV.
*
* Returns the inverse of the left-tailed probability of the chi-squared distribution.
*
* @param mixed (float) $probability Probability for the function
* @param mixed (int) $degrees degrees of freedom
*
* @return float|string
*/
public static function inverseLeftTail($probability, $degrees)
{
$probability = Functions::flattenSingleValue($probability);
$degrees = Functions::flattenSingleValue($degrees);

try {
$probability = self::validateFloat($probability);
$degrees = self::validateInt($degrees);
} catch (Exception $e) {
return $e->getMessage();
}

if ($probability < 0.0 || $probability > 1.0 || $degrees < 1) {
return Functions::NAN();
}

return self::inverseLeftTailCalculation($probability, $degrees);
}

public static function test($actual, $expected)
{
$rows = count($actual);
Expand Down Expand Up @@ -167,4 +198,104 @@ protected static function degrees(int $rows, int $columns): int

return ($columns - 1) * ($rows - 1);
}

private static function inverseLeftTailCalculation($probability, $degrees)
{
// bracket the root
$min = 0;
$sd = sqrt(2.0 * $degrees);
$max = 2 * $sd;
$s = -1;

while ($s * self::pchisq($max, $degrees) > $probability * $s) {
$min = $max;
$max += 2 * $sd;
}

// Find root using bisection
$chi2 = 0.5 * ($min + $max);

while (($max - $min) > self::EPS * $chi2) {
if ($s * self::pchisq($chi2, $degrees) > $probability * $s) {
$min = $chi2;
} else {
$max = $chi2;
}
$chi2 = 0.5 * ($min + $max);
}

return $chi2;
}

private static function pchisq($chi2, $degrees)
{
return self::gammp($degrees, 0.5 * $chi2);
}

private static function gammp($n, $x)
{
if ($x < 0.5 * $n + 1) {
return self::gser($n, $x);
}

return 1 - self::gcf($n, $x);
}

// Return the incomplete gamma function P(n/2,x) evaluated by
// series representation. Algorithm from numerical recipe.
// Assume that n is a positive integer and x>0, won't check arguments.
// Relative error controlled by the eps parameter
private static function gser($n, $x)
{
$gln = Gamma::ln($n / 2);
$a = 0.5 * $n;
$ap = $a;
$sum = 1.0 / $a;
$del = $sum;
for ($i = 1; $i < 101; ++$i) {
++$ap;
$del = $del * $x / $ap;
$sum += $del;
if ($del < $sum * self::EPS) {
break;
}
}

return $sum * exp(-$x + $a * log($x) - $gln);
}

// Return the incomplete gamma function Q(n/2,x) evaluated by
// its continued fraction representation. Algorithm from numerical recipe.
// Assume that n is a postive integer and x>0, won't check arguments.
// Relative error controlled by the eps parameter
private static function gcf($n, $x)
{
$gln = Gamma::ln($n / 2);
$a = 0.5 * $n;
$b = $x + 1 - $a;
$fpmin = 1.e-300;
$c = 1 / $fpmin;
$d = 1 / $b;
$h = $d;
for ($i = 1; $i < 101; ++$i) {
$an = -$i * ($i - $a);
$b += 2;
$d = $an * $d + $b;
if (abs($d) < $fpmin) {
$d = $fpmin;
}
$c = $b + $an / $c;
if (abs($c) < $fpmin) {
$c = $fpmin;
}
$d = 1 / $d;
$del = $d * $c;
$h = $h * $del;
if (abs($del - 1) < self::EPS) {
break;
}
}

return $h * exp(-$x + $a * log($x) - $gln);
}
}
Original file line number Diff line number Diff line change
@@ -0,0 +1,37 @@
<?php

namespace PhpOffice\PhpSpreadsheetTests\Calculation\Functions\Statistical;

use PhpOffice\PhpSpreadsheet\Calculation\Functions;
use PhpOffice\PhpSpreadsheet\Calculation\Statistical;
use PHPUnit\Framework\TestCase;

class ChiInvLeftTailTest extends TestCase
{
protected function setUp(): void
{
Functions::setCompatibilityMode(Functions::COMPATIBILITY_EXCEL);
}

/**
* @dataProvider providerCHIINV
*
* @param mixed $expectedResult
* @param mixed $probability
* @param mixed $degrees
*/
public function testCHIINV($expectedResult, $probability, $degrees): void
{
$result = Statistical\Distributions\ChiSquared::inverseLeftTail($probability, $degrees);
if (!is_string($expectedResult)) {
$reverse = Statistical\Distributions\ChiSquared::distributionLeftTail($result, $degrees, true);
self::assertEqualsWithDelta($probability, $reverse, 1E-12);
}
self::assertEqualsWithDelta($expectedResult, $result, 1E-12);
}

public function providerCHIINV()
{
return require 'tests/data/Calculation/Statistical/CHIINVLeftTail.php';
}
}
64 changes: 64 additions & 0 deletions tests/data/Calculation/Statistical/CHIINVLeftTail.php
Original file line number Diff line number Diff line change
@@ -0,0 +1,64 @@
<?php

return [
[
4.671330448981,
0.3, 7,
],
[
12.548861396889,
0.75, 10,
],
[
3.283020286760,
0.93, 1,
],
[
1.832581463748,
0.6, 2,
],
[
0.454936423120,
0.5, 1,
],
[
4.351460191096,
0.5, 5,
],
[
1.323303696932,
0.75, 1,
],
[
0.210721031316,
0.1, 2,
],
[
3.218875824869,
0.8, 2,
],
[
1.212532903046,
0.25, 3,
],
[
'#VALUE!',
'NaN', 3,
],
[
'#VALUE!',
0.25, 'NaN',
],
'Probability < 0' => [
'#NUM!',
-0.1, 3,
],
'Probability > 1' => [
'#NUM!',
1.1, 3,
],
'Freedom > 1' => [
'#NUM!',
0.1, 0.5,
],
];

0 comments on commit e68978f

Please sign in to comment.