寻求Delphi版Chi-square分布可用代码及浮点溢出问题解决

Let's break down why your code hits overflow with large chi-square values, and how to fix it without relying on third-party libraries.

The Root Cause #

Your Reihe function calculates each summand directly as power(chi, 2*k) / Bruch, where Bruch is the product of (f + 2*i) for i=1 to k. When chi is large (like 138.6), chi^(2*k) grows exponentially far faster than the denominator, quickly exceeding the maximum value representable by Delphi's Real type—this triggers the floating-point overflow.

The Solution: Recursive Summand Calculation #

Instead of computing each summand from scratch (which leads to huge intermediate values), we can calculate each term recursively using the previous term. This keeps numbers manageable and avoids overflow. Here's the recurrence logic:

• The first summand (k=1) is chi² / (f + 2)
• Each subsequent summand is previous_summand * chi² / (f + 2*k)

We also added a relative error check to stop the loop early if the summand becomes negligible compared to the total sum, which is more robust for large datasets than an absolute threshold.

Modified Full Code #

function LnGamma(x: Real): Real;
const
  a0 = 0.083333333096;
  a1 = -0.002777655457;
  a2 = 0.000777830670;
  c = 0.918938533205;
var
  r: Real;
begin
  r := (a0 + (a1 + a2 / Sqr(x)) / Sqr(x)) / x;
  LnGamma := (x - 0.5) * Ln(x) - x + c + r;
end;

function LnFak(x: Real): Real;
var
  z: Real;
begin
  z := x + 1;
  LnFak := LnGamma(z);
end;

function Reihe(chi: Real; f: Real): Real;
const
  MaxError = 0.0001;
var
  Summe, Summand, chiSquared: Real;
  k: Integer;
begin
  Summe := 1.0;
  if chi = 0 then
  begin
    Result := Summe;
    Exit;
  end;

  chiSquared := Sqr(chi);
  // Initialize first summand (k=1)
  Summand := chiSquared / (f + 2);
  Summe := Summe + Summand;
  k := 2;

  repeat
    // Calculate next summand using previous term
    Summand := Summand * chiSquared / (f + 2*k);
    Summe := Summe + Summand;
    Inc(k);
    // Stop when summand is negligible relative to the sum
  until (Abs(Summand) < MaxError * Abs(Summe)) or (Summand = 0);

  Result := Summe;
end;

function IntegralChi(chisqr: Real; f: Integer): Real;
var
  s, halfChisqr, lnTerm: Real;
begin
  if chisqr <= 0 then
  begin
    Result := 1.0;
    Exit;
  end;

  halfChisqr := 0.5 * chisqr;
  // Calculate the term using logarithms to avoid overflow in intermediate steps
  lnTerm := (f/2) * Ln(halfChisqr) - halfChisqr - LnGamma((f + 2)/2);
  s := Exp(lnTerm) * Reihe(Sqrt(chisqr), f);
  Result := 1 - s;
end;

Key Improvements #

1. Recursive Summand Calculation: The Reihe function now builds each term from the previous one, preventing exponential growth of intermediate values.
2. Relative Error Check: The loop stops when the summand is smaller than MaxError * Summe, which is more reliable for large sums than an absolute threshold.
3. Logarithm-Based Term Calculation: In IntegralChi, we compute the exponential term using logarithms first, then exponentiate once—this avoids overflow from multiplying large numbers directly.
4. Edge Case Handling: Added checks for chi=0 and chisqr<=0 to return correct values immediately.

Testing the Fix #

For your problematic case (Chi=138.609137, df=4), the modified code returns a p-value extremely close to 0 (effectively 0 for practical purposes) without hitting a floating-point overflow. For your working case (Chi=1.142132, df=1), it still returns ~0.285202, matching SPSS results.

内容的提问来源于stack exchange，提问作者kwadratens