代码之家 › 专栏 › 技术社区 › kwadratens

Delphi中的卡方分布函数代码

chi-squared distribution statistics delphi

1

kwadratens · 技术社区 · 6 年前

我一直在寻找可用的完整代码 chi-square 德尔福的分销。有些代码是通过网络编写的,但通常不起作用或缺少部分,无法编译等。。还有一些库,但我对一些可以简单实现的代码感兴趣。

我发现有些东西几乎能起作用。一些德语部分已经修复,它编译并给出 p-values 对于大多数数据:

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;

它对相对较大的结果非常有效。

例如:

对于 Chi = 1.142132 和 df = 1 我要走了 p 关于 0.285202 ,这是完美的。等同于 SPSS 结果或其他程序。

但例如 Chi = 138.609137 和 df = 4 我应该收到一些关于 0.000000 ,但中出现浮点溢出错误 Reiche 作用 Summe 和 Summand 那是非常大的。

我承认理解分布函数不是我的强项,所以也许有人会告诉我我做错了什么?

非常感谢您提供的信息

1 回复 | 直到 6 年前

1

4

gammatester 6 年前

您应该调试您的程序并发现有溢出在k=149的循环中。k=148时,Bruch值为3.3976725289e+304。下一步计算Bruch溢出。解决方法是编写代码

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

有了这个变化,你就得到了价值 IntegralChi(138.609137,4) = 1.76835197E-7 经过156次迭代。

请注意,您的计算(即使对于这个简单的算法)是次优的因为你反复计算布鲁赫值。只需更新一次每个循环:

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;

类似的考虑也应适用于计算 power(chi, 2*k) 然后将其与改进的Bruch评估相结合。

编辑: 作为对您评论的回应,下面是基于幂函数特性的改进版本,即 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;