Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Chi squared inverse left tailed #1964

Merged
merged 12 commits into from
Mar 28, 2021
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,
],
];