CODEWITHSUNDEEP
delphistatisticsdistributionchi-squared
kwadratens
kwadratens

Reputation: 345

Code for Chi-square distribution function in Delphi

I have been looking for usable and full code for chi-square distribution in Delphi. There are some codes via net, but usually they don't work or have missing parts, do not compile etc.. There are also some libraries, but I'm interested about some code that I just can simply implement.

I've found something almost working. Some german parts have been fixed, it compiles and it gives p-values for most of the data:

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
  Bruch,
  Summe,
  Summand : Real;
  k, i    : longint;    
begin
  Summe := 1;
  k := 1;
  repeat
    Bruch := 1;
    for i := 1 to k do
      Bruch := Bruch * (f + 2 * i);
    Summand := power(chi, 2 * k) / Bruch;
    Summe := Summe + Summand;
    k := succ(k);
  until (Summand < MaxError);
  Reihe := Summe;
end;

function IntegralChi (chisqr : Real; f : longint) : Real;
var
  s : Real;
begin
  S := power((0.5 * chisqr), f/2) * Reihe(sqrt(chisqr), f)
                  * exp((-chisqr/2) - LnGamma((f + 2) / 2));
  IntegralChi := 1 - s;
end;

It works quite good for relatively big results.

For example:

For Chi = 1.142132 and df = 1 I'm getting p about 0.285202, which is perfect. Same as SPSS result or other programs.

But for example Chi = 138.609137 and df = 4 I should recieive something about 0.000000, but I'm getting floating point overflow error in Reiche function. Summe and Summand are very big then.

I admit that understanding distribution function is not my strong point, so maybe someone will tell me what I did wrong?

Thank you very much for the information

Upvotes: 1

Views: 350

Answers (1)

gammatester
gammatester

Reputation: 1141

You should debug your program and find that there is an overflow in your loop for k=149. For k=148 the value of Bruch is 3.3976725289e+304. The next computation of Bruch overflows. A fix is to code

for i := 1 to k do
  Bruch := Bruch / (f + 2 * i);
Summand := power(chi, 2 * k) * Bruch;

With this change you get the value IntegralChi(138.609137,4) = 1.76835197E-7 after 156th iteration.

Note that your computation (even for this simple algorithm) is sub-optimal because you compute the Bruch value over and over again. Just update it once per loop:

function Reihe (chi : Real; f : Real) : Real;
  const MaxError = 0.0001;
var
  Bruch,
  Summe,
  Summand : Real;
  k    : longint;
begin
  Summe := 1;
  k := 1;
  Bruch := 1;
  repeat
    Bruch := Bruch / (f + 2 * k);
    Summand := power(chi, 2 * k) * Bruch;
    Summe := Summe + Summand;
    k := succ(k);
  until (Summand < MaxError);
  Reihe := Summe;
end;

Similar consideration should be applied to compute power(chi, 2*k) and then combine this with the improved evaluation of Bruch.

Edit: As a response to your comment, here the improved version based on the property of the power function, that is power(chi, 2*(k+1)) = power(chi, 2*k)*sqr(chi)

function Reihe (chi : Real; f : Real) : Real;
  const MaxError = 0.0001;
var
  chi2,
  Summe,
  Summand : Real;
  k    : longint;
begin
  Summe := 1;
  k := 1;
  Summand := 1;
  chi2 := sqr(chi);
  repeat
    Summand := Summand * chi2 / (f + 2 * k);
    Summe := Summe + Summand;
    k := succ(k);
  until (Summand < MaxError);
  Reihe := Summe;
end;

Upvotes: 4

Related Questions