geckon01/PhpSpreadsheet
1
0
Fork 0
mirror of https://github.com/PHPOffice/PhpSpreadsheet.git synced 2026-08-30 03:59:11 +00:00
Code Issues Packages Projects Releases Wiki Activity

Fix incomplete gamma convergence for GAMMA.DIST / CHISQ.DIST

Browse Source 
The incomplete gamma primitive used a fixed 32-term power series with no
convergence test, so GAMMA.DIST, GAMMADIST, CHISQ.DIST(.RT), GAMMAINV and
CHISQ.INV were grossly wrong once the series argument reached ~32
(e.g. CHISQ.DIST.RT(80, 4) returned 0.806 instead of 1.74e-16).

Replace it with the standard convergence-tested regularized incomplete
gamma: series P(a,x) for x < a+1, continued fraction Q(a,x) for x >= a+1.
CHISQ.DIST.RT now uses Q directly so the right tail stays free of
1 - P cancellation. Consolidates the duplicate copy that already existed
privately in ChiSquared onto the shared primitive.
This commit is contained in:
gaoflow
2026-07-28 09:10:13 +02:00
committed by Vincent Gao
parent 8b4b934011
commit e88cd615cc
8 changed files with 134 additions and 88 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/PhpSpreadsheet/Calculation/Statistical/Distributions/ChiSquared.php
+3 -73
View File
@@ -50,7 +50,7 @@ class ChiSquared
return ExcelError::NAN();
}
return 1 - (Gamma::incompleteGamma($degrees / 2, $value / 2) / Gamma::gammaValue($degrees / 2));
return Gamma::regularizedGammaQ($degrees / 2, $value / 2);
}
/**
@@ -133,8 +133,7 @@ class ChiSquared
return ExcelError::NAN();
}
$callback = fn (float $value): float => 1 - (Gamma::incompleteGamma($degrees / 2, $value / 2)
/ Gamma::gammaValue($degrees / 2));
$callback = fn (float $value): float => Gamma::regularizedGammaQ($degrees / 2, $value / 2);
$newtonRaphson = new NewtonRaphson($callback);
@@ -258,75 +257,6 @@ class ChiSquared
private static function pchisq(float $chi2, int $degrees): float
{
return self::gammp($degrees, 0.5 * $chi2);
}
private static function gammp(int $n, float $x): float
{
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(int $n, float $x): float
{
/** @var float $gln */
$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(int $n, float $x): float
{
/** @var float $gln */
$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);
return Gamma::regularizedGammaP($degrees / 2, 0.5 * $chi2);
}
}
src/PhpSpreadsheet/Calculation/Statistical/Distributions/GammaBase.php
+83 -10
View File
@@ -20,7 +20,7 @@ abstract class GammaBase
protected static function calculateDistribution(float $value, float $a, float $b, bool $cumulative): float
{
if ($cumulative) {
return self::incompleteGamma($a, $value / $b) / self::gammaValue($a);
return self::regularizedGammaP($a, $value / $b);
}
return (1 / ($b ** $a * self::gammaValue($a))) * $value ** ($a - 1) * exp(0 - ($value / $b));
@@ -96,17 +96,90 @@ abstract class GammaBase
//
public static function incompleteGamma(float $a, float $x): float
{
static $max = 32;
$summer = 0;
for ($n = 0; $n <= $max; ++$n) {
$divisor = $a;
for ($i = 1; $i <= $n; ++$i) {
$divisor *= ($a + $i);
}
$summer += ($x ** $n / $divisor);
// Unregularized lower incomplete gamma; kept for backward compatibility.
return self::regularizedGammaP($a, $x) * self::gammaValue($a);
}
/**
* Regularized lower incomplete gamma P(a,x) = gamma(a,x) / Gamma(a).
* Series for x < a+1, else the complement of the continued fraction.
*/
public static function regularizedGammaP(float $a, float $x): float
{
if ($x <= 0.0 || $a <= 0.0) {
return 0.0;
}
if ($x < $a + 1.0) {
return self::gammaSeries($a, $x);
}
return $x ** $a * exp(0 - $x) * $summer;
return 1.0 - self::gammaContinuedFraction($a, $x);
}
/**
* Regularized upper incomplete gamma Q(a,x) = 1 - P(a,x).
* Continued fraction for x >= a+1 keeps the right tail free of cancellation.
*/
public static function regularizedGammaQ(float $a, float $x): float
{
if ($x <= 0.0 || $a <= 0.0) {
return 1.0;
}
if ($x < $a + 1.0) {
return 1.0 - self::gammaSeries($a, $x);
}
return self::gammaContinuedFraction($a, $x);
}
// P(a,x) by its series representation (Numerical Recipes gser).
private static function gammaSeries(float $a, float $x): float
{
$gln = self::logGamma($a);
$ap = $a;
$sum = 1.0 / $a;
$del = $sum;
for ($i = 1; $i <= self::MAX_ITERATIONS; ++$i) {
++$ap;
$del *= $x / $ap;
$sum += $del;
if (abs($del) < abs($sum) * self::EPS) {
break;
}
}
return $sum * exp(-$x + $a * log($x) - $gln);
}
// Q(a,x) by its continued fraction representation (Numerical Recipes gcf).
private static function gammaContinuedFraction(float $a, float $x): float
{
$fpMin = 1.0e-300;
$gln = self::logGamma($a);
$b = $x + 1.0 - $a;
$c = 1.0 / $fpMin;
$d = 1.0 / $b;
$h = $d;
for ($i = 1; $i <= self::MAX_ITERATIONS; ++$i) {
$an = -$i * ($i - $a);
$b += 2.0;
$d = $an * $d + $b;
if (abs($d) < $fpMin) {
$d = $fpMin;
}
$c = $b + $an / $c;
if (abs($c) < $fpMin) {
$c = $fpMin;
}
$d = 1.0 / $d;
$del = $d * $c;
$h *= $del;
if (abs($del - 1.0) < self::EPS) {
break;
}
}
return $h * exp(-$x + $a * log($x) - $gln);
}
private const GAMMA_VALUE_P0 = 1.000000000190015;
tests/PhpSpreadsheetTests/Calculation/Functions/Statistical/ChiDistRightTailTest.php
+21
View File
@@ -19,6 +19,27 @@ class ChiDistRightTailTest extends AllSetupTeardown
return require 'tests/data/Calculation/Statistical/CHIDISTRightTail.php';
}
// Deep right tail computed directly via Q(a,x); "1 - P" would collapse these to 0.
// Expected values from mpmath (gammainc regularized, dps 30 == 35).
#[\PHPUnit\Framework\Attributes\DataProvider('providerChiDistRightTailDeepTail')]
public function testChiDistRightTailDeepTail(float $expectedResult, float $value, int $degrees): void
{
$calculation = Calculation::getInstance();
$formula = "=CHISQ.DIST.RT($value, $degrees)";
/** @var float $result */
$result = $calculation->calculateFormula($formula);
self::assertEqualsWithDelta($expectedResult, $result, abs($expectedResult) * 1.0e-9, $formula);
}
public static function providerChiDistRightTailDeepTail(): array
{
return [
[1.7418252446695515e-16, 80.0, 4],
[5.0600460658425739e-21, 120.0, 10],
[1.1253473960842734e-31, 200.0, 20],
];
}
#[\PHPUnit\Framework\Attributes\DataProvider('providerChiDistRightTailArray')]
public function testChiDistRightTailArray(array $expectedResult, string $values, string $degrees): void
{
tests/PhpSpreadsheetTests/Calculation/Functions/Statistical/ChiInvRightTailTest.php
+2 -2
View File
@@ -49,8 +49,8 @@ class ChiInvRightTailTest extends AllSetupTeardown
return [
'row/column vectors' => [
[
[7.8061229155968075, 6.345811195521517, 100.0],
[13.266097125199911, 11.34032237742413, 24.388802783239434],
[7.8061229155968075, 6.345811195521517, 16.907871682617596],
[13.266097125199924, 11.340322377424133, 24.388802639997067],
],
'{0.35, 0.5, 0.018}',
'{7; 12}',
tests/PhpSpreadsheetTests/Calculation/Functions/Statistical/GammaInvTest.php
+2 -2
View File
@@ -76,8 +76,8 @@ class GammaInvTest extends AllSetupTeardown
'row/column vectors' => [
[
[2.772588722239782, 5.38526905777939, 12.548861396889375],
[5.545177444479563, 10.77053811555878, 25.09772279377875],
[6.931471805599453, 13.463172644448473, 31.372153492223436],
[5.545177444479561, 10.770538115558788, 25.097722793778765],
[6.931471805599456, 13.463172644448484, 31.372153492223447],
],
'0.75',
'{1, 2, 5}',
tests/data/Calculation/Statistical/CHIDISTRightTail.php
+9
View File
@@ -35,6 +35,15 @@ return [
0.046011705689,
8, 3,
],
// Large x: the old fixed 32-term series diverged badly here.
[
0.177045452100,
70, 60,
],
[
0.064570368921,
100, 80,
],
[
'#VALUE!',
'NaN', 3,
tests/data/Calculation/Statistical/CHITEST.php
+1 -1
View File
@@ -9,7 +9,7 @@ return [
[[100, 85], [70, 90]],
],
[
0.000308192027,
0.0003081920170083,
[[58, 35], [11, 25], [10, 23]],
[[45.35, 47.65], [17.56, 18.44], [16.09, 16.91]],
],
tests/data/Calculation/Statistical/GAMMADIST.php
+13
View File
@@ -23,6 +23,19 @@ return [
0.576809918873,
6, 3, 2, 1,
],
// Large x/b: the old fixed 32-term series diverged badly here.
'Large x mid-range' => [
0.522481190430,
35, 35, 1, true,
],
'Large x saturating' => [
1.0,
50, 2, 1, true,
],
[
1.0,
80, 2, 1, true,
],
[
'#VALUE!',
'NAN', 3, 2, true,