Функция распределения Пирсона

← →
Jimmy (2009-10-26 19:07) [0]

Считаю значение функции распределения Пирсона в точке как сумму степенного ряда. Проблема в том, что среди членов ряда есть очень большие числа, порядка 10^30, но сумма меньше 1. Происходит потеря точности в каждом слагаемом что приводит к совершенно нелепому результату. Подскажите пожалуйста как можно решить проблему?

← →
Юрий Зотов © (2009-10-26 19:10) [1]

Группируйте члены ряда так, чтобы каждая частичная сумма была близка к 1, а потом складывайте эти частичные суммы.

← →
Юрий Зотов © (2009-10-26 19:12) [2]

Пример (условный):

10^30 + 1 - 10^30 - получим ноль
10^30 - 10^30 + 1 - получим единицу

← →
palva © (2009-10-26 22:04) [3]

> как можно решить проблему?

Если проблема просто посчитать, то лучше использовать готовый алгоритм - найти в какой-нибудь книжке. Обычно такой алгоритм использует аппроксимацию в тех местах, где ряд плохо сходится.

Вот здесь поищите http://www.srcc.msu.su/num_anal/lib_na/cat/cat66.htm
Программы на фортране и си.

Распределение Пирсона ( оно же хи-квадрат) обычно считается через неполный интеграл гамма-функции.
См. ниже:

const
MathEps = 1.1102230246251565E-16; // 2**-53}
MaxLog = 7.0978271289338399E2; // log(2**1024)}
MinLog = -7.08396418532264106224E2; // log(2**-1022}
Ln2 = 6.9314718055994531E-1; // log(2)}
pio4 = 0.7853981633974483E0; // Pi/4}
tpi = 0.6366197723675813E0; // 2/Pi
Sqrt2 = 1.41421356237309504880; // sqrt(2)}
Infinity = MaxDouble;

function ChDtr(df: Double; x: Double): Double; // Хи-квадрат распределение
//ChDtr( v, x ) = IGam( v/2.0, x/2.0 )

function IGam(a, x: Double): Double; // Неполный интеграл Гамма функции

function IGamC(a, x: Double): Double; // Дополнительный неполный интеграл Гамма функции
// 0<x<100, 0<a<100

function LnGam(X: Double): Double; // Натуральный логарифм Гамма функции
// |arg|<2.55e305

implementation

procedure MathError(TypeError: TTypeError;
NameFunc: string);
var S: string;
begin
case TypeError of
Domain: S := "Аргумент " + NameFunc + " вне области определения";
Sing: S := "Функция " + NameFunc + " в точке сингулярности (разрыв)";
Overflow: S := NameFunc + " - переполнение (overflow)";
Underflow: S := NameFunc + " - антипереполнение (underflow) ";
TLoss, PLoss: S := NameFunc + " - потеря точности вычислений";
end;
raise EFuncLib.Create(S);
end;

function chdtr(df: Double; x: Double): Double;
begin
if (x < 0.0) or (df < 1.0) then begin MathError(DOMAIN, "chdtr"); Result := 0; Exit; end;
Result := Igam(df / 2.0, x / 2.0);
end;

function IGam(a, x: Double): Double;
var ans, ax, c, r: Double;
begin
if (x <= 0) or (a <= 0) then begin igam := 0.0; Exit; end;
if (x > 1.0) and (x > a) then begin igam := (1.0 - igamc(a, x)); Exit; end;
ax := (a * Ln(x)) - x - LnGam(a);
if ax < -MAXLOG then begin MathError(UNDERFLOW, "igam"); igam := 0; Exit; end;
ax := Exp(ax);
r := a; c := 1.0; ans := 1.0;
repeat
r := r + 1.0; c := c * x / r; ans := ans + c;
until (c / ans <= MathEps);
IGam := (ans * ax / a);
end;

function IGamC(a, x: Double): Double;
var ans, ax, c, yc, r, t, y, z, pk, pkm1, pkm2, qk, qkm1, qkm2: Double;
begin
if (x <= 0) or (a <= 0) then begin igamc := 1.0; Exit; end;
if (x < 1.0) or (x < a) then begin igamc := (1.0 - igam(a, x)); Exit; end;
ax := (a * Ln(x)) - x - LnGam(a);
if (ax < -MAXLOG) then begin MathError(UNDERFLOW, "igamc"); igamc := 0.0; exit end;
ax := Exp(ax);
{ continued fraction }
y := 1.0 - a; z := x + y + 1.0; c := 0.0; pkm2 := 1.0; qkm2 := x; pkm1 := x + 1.0;
qkm1 := z * x; ans := pkm1 / qkm1;
repeat
c := c + 1.0; y := y + 1.0; z := z + 2.0; yc := y * c; pk := (pkm1 * z) - (pkm2 * yc);
qk := (qkm1 * z) - (qkm2 * yc);
if qk <> 0 then begin r := pk / qk; t := Abs((ans - r) / r); ans := r; end
else t := 1.0;
pkm2 := pkm1; pkm1 := pk; qkm2 := qkm1; qkm1 := qk;
if Abs(pk) > big then
begin
pkm2 := pkm2 * biginv; pkm1 := pkm1 * biginv;
qkm2 := qkm2 * biginv; qkm1 := qkm1 * biginv;
end;
until t <= MathEps;
Result := (ans * ax);
end;

function LnGam(x: Double): Double;
const
LOGPI = 1.14472988584940017414;
LS2PI = 0.91893853320467274178;
MAXLGM = 2.556348E305;
A: array[0..4] of Double =
(8.11614167470508450300E-4, -5.95061904284301438324E-4,
7.93650340457716943945E-4, -2.77777777730099687205E-3,
8.33333333333331927722E-2);
B: array[0..5] of Double =
(-1.37825152569120859100E3, -3.88016315134637840924E4,
-3.31612992738871184744E5, -1.16237097492762307383E6,
-1.72173700820839662146E6, -8.53555664245765465627E5);
C: array[0..5] of Double =
(-3.51815701436523470549E2, -1.70642106651881159223E4,
-2.20528590553854454839E5, -1.13933444367982507207E6,
-2.53252307177582951285E6, -2.01889141433532773231E6);
var p, q, u, w, z: Double;
i: Integer;
begin
SgnGam := 1;
if x < -34.0 then
begin
q := -x;
w := LnGam(q); {note this modifies sgngam!}
p := floor(q);
if p = q then
begin MathError(Sing, "LnGam"); Result := INFINITY; Exit; end;
i := Round(p);
if (i and 1) = 0 then sgngam := -1 else sgngam := 1;
z := q - p;
if (z > 0.5) then begin p := p + 1.0; z := p - q; end;
z := q * Sin(Pi * z);
if z = 0.0 then begin MathError(Sing, "LnGam"); Result := INFINITY; Exit; end;
Result := LOGPI - Ln(z) - w;
end;
if x < 13.0 then
begin
z := 1.0; p := 0.0; u := x;
while u >= 3.0 do begin p := p - 1.0; u := x + p; z := z * u; end;
while (u < 2.0) do
begin
if u = 0.0 then begin MathError(Sing, "LnGam"); Result := INFINITY; Exit; end;
z := z / u; p := p + 1.0; u := x + p;
end;
if z < 0.0 then begin sgngam := -1; z := -z; end
else sgngam := 1;
if (u = 2.0) then begin Result := (Ln(z)); Exit; end;
p := p - 2.0; x := x + p;
p := x * PolEvl(x, B) / p1evl(x, C);
Result := (Ln(z) + p);
Exit;
end { if };
if (x > 2.556348E305) then begin Result := SgnGam * INFINITY; Exit; end;
q := ((x - 0.5) * Ln(x)) - x + LS2PI;
if x > 1.0E8 then begin Result := q; Exit; end;
p := 1.0 / (x * x);
if (x >= 1000.0) then
q := q + ((((7.9365079365079365079365E-4 * p) -
2.7777777777777777777778E-3) * p) +
0.0833333333333333333333) / x
else q := q + PolEvl(p, A) / x;
Result := (q);
end;

function PolEvl(X: Extended; const Coefficients: array of Double): Extended;
{ Horner"s method }
var
I: Integer;
begin
Result := Coefficients[Low(Coefficients)];
for I := Low(Coefficients) + 1 to High(Coefficients) do
Result := Result * X + Coefficients[I];
end;

Функция распределения Пирсона Найти похожие ветки