1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
|
{ ******************************************************************
Inverses of incomplete Gamma function and Khi-2 distribution
Translated from C code in Cephes library (http://www.moshier.net)
****************************************************************** }
unit uinvgam;
interface
uses
utypes, ugamma, uigamma, uinvnorm;
function InvGamma(A, P : Float) : Float;
{ ------------------------------------------------------------------
Given P, the function finds X such that
IGamma(A, X) = P
It is best valid in the right-hand tail of the distribution, P > 0.5
------------------------------------------------------------------ }
function InvKhi2(Nu : Integer; P : Float) : Float;
{ ------------------------------------------------------------------
Inverse of Khi-2 distribution function
Returns the argument, X, for which the area under the Khi-2
probability density function (integrated from 0 to X)
is equal to P.
------------------------------------------------------------------ }
implementation
function InvGamma(A, P : Float) : Float;
var
x0, x1, x, Y, yl, yh, y1, d, lgm, dithresh : Float;
i, ndir : Integer;
label
ihalve, cont, cont1, done;
begin
if P > 0.5 then SetErrCode(FPLoss) else SetErrCode(FOk);
Y := 1.0 - P;
{ Bound the solution }
x0 := MaxNum;
yl := 0.0;
x1 := 0.0;
yh := 1.0;
dithresh := 5 * MachEp;
{ Approximation to inverse function }
d := 1.0 / (9.0 * a);
y1 := 1.0 - d - InvNorm(Y) * sqrt(d);
x := a * y1 * y1 * y1;
lgm := LnGamma(a);
for i := 0 to 9 do
begin
if (x > x0) or (x < x1) then goto ihalve;
y1 := JGamma(a, x);
if (y1 < yl) or (y1 > yh) then goto ihalve;
if y1 < Y then
begin
x0 := x;
yl := y1
end
else
begin
x1 := x;
yh := y1
end;
{ Compute the derivative of the function at this point }
d := (a - 1) * Ln(x) - x - lgm;
if d < MinLog then goto ihalve;
d := -exp(d);
{ Compute the step to the next approximation of x }
d := (y1 - Y) / d;
if abs(d / x) < MachEp then goto done;
x := x - d;
end;
{ Resort to interval halving if Newton iteration did not converge }
ihalve:
d := 0.0625;
if x0 = MaxNum then
begin
if x <= 0 then x := 1;
while x0 = MaxNum do
begin
x := (1 + d) * x;
y1 := JGamma(a, x);
if y1 < Y then
begin
x0 := x;
yl := y1;
goto cont
end;
d := d + d
end
end;
cont:
d := 0.5;
ndir := 0;
for i := 0 to 399 do
begin
x := x1 + d * (x0 - x1);
y1 := JGamma(a, x);
lgm := (x0 - x1) / (x1 + x0);
if abs(lgm) < dithresh then goto cont1;
lgm := (y1 - Y) / Y;
if abs(lgm) < dithresh then goto cont1;
if x <= 0 then goto cont1;
if y1 >= Y then
begin
x1 := x;
yh := y1;
if ndir < 0 then
begin
ndir := 0;
d := 0.5
end
else if ndir > 1 then
d := 0.5 * d + 0.5
else
d := (Y - yl) / (yh - yl);
ndir := ndir + 1;
end
else
begin
x0 := x;
yl := y1;
if ndir > 0 then
begin
ndir := 0;
d := 0.5
end
else if ndir < -1 then
d := 0.5 * d
else
d := (Y - yl) / (yh - yl);
ndir := ndir - 1;
end;
end;
cont1:
if x = 0 then SetErrCode(FUnderflow);
done:
InvGamma := x
end;
function InvKhi2(Nu : Integer; P : Float) : Float;
{ ------------------------------------------------------------------
Inverse of Khi-2 distribution function
Returns the argument, X, for which the area under the Khi-2
probability density function (integrated from 0 to X)
is equal to P.
------------------------------------------------------------------ }
begin
if (P < 0.0) or (P > 1.0) or (Nu < 1) then
InvKhi2 := DefaultVal(FDomain, 0.0)
else
InvKhi2 := 2.0 * InvGamma(0.5 * Nu, P);
end;
end.
|
