Код для функции распределения хи-квадрат в Delphi - PullRequest
Новая
       59

Код для функции распределения хи-квадрат в Delphi

0 голосов
/ 17 мая 2018

Я искал полезный и полный код для распространения chi-square в Delphi.Есть некоторые коды через сеть, но обычно они не работают или имеют недостающие части, не компилируются и т. Д. Есть также некоторые библиотеки, но меня интересует некоторый код, который я просто могу реализовать.

Я нашел что-то почти работающее.Некоторые немецкие части были исправлены, он компилируется и выдает 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 Ответ

0 голосов
/ 17 мая 2018

Вы должны отладить свою программу и обнаружить переполнение в вашем цикле для к = 149. Для k = 148 значение Бруха составляет 3.3976725289e + 304. Следующее вычисление переполнения Бруха. Исправление в коде 

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), а затем объединить его с улучшенной оценкой Бруха.

Редактировать: В качестве ответа на ваш комментарий, здесь улучшенная версия, основанная на свойстве степенной функции, то есть 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;
...