cdfchii.src
来自「没有说明」· SRC 代码 · 共 563 行 · 第 1/2 页
SRC
563 行
/*
** cdfchii.src
**
**
** (C) Copyright 1995-1996 Aptech Systems, Inc.
** All Rights Reserved.
**
** This Software Product is PROPRIETARY SOURCE CODE OF APTECH
** SYSTEMS, INC. This File Header must accompany all files using
** any portion, in whole or in part, of this Source Code. In
** addition, the right to create such files is strictly limited by
** Section 2.A. of the GAUSS Applications License Agreement
** accompanying this Software Product.
**
** If you wish to distribute any portion of the proprietary Source
** Code, in whole or in part, you must first obtain written
** permission from Aptech Systems.
**
**-------------------**------------------**-------------------**-----------**
**-------------------**------------------**-------------------**-----------**
**
**
** Format Purpose Line
** ---------------------------------------------------------------------------
** c = cdfchii(p,n); inverse chi-square 30
** x = gammaii(a,p) inverse incomplete gamma 59
*/
/*
**> cdfchii
**
** Purpose: compute chi-square abscissae values given probability
** and degrees of freedom.
**
** Format: c = cdfchii(p,n);
**
** Input: p MxN matrix, or Mx1, or 1xN vector, probabilities.
**
** n MxN matrix, or Mx1, or 1xN vector, degrees of freedom.
**
** Output: c MxN matrix, abscissae values for chi-square
** distribution.
**
** Globals: gammaii
*/
proc cdfchii(p,n);
if p < 0 or p > 1 or n <= 0;
if not trapchk(1);
errorlog "CDFCHII: argument(s) out of range";
end;
endif;
retp(error(0));
endif;
retp(2*gammaii(0.5*n,p));
endp;
/*
**> gammaii
**
** Purpose: Compute the inverse incomplete Gamma function.
**
** Format: x = gammaii(a,p);
**
** Input: a MxN matrix, or Mx1, or 1xN vector, exponents.
**
** p MxN matrix, or Mx1, or 1xN vector, incomplete gamma
** values.
**
** Output: x MxN matrix, abscissae.
**
** Globals: _ginvinc, __macheps
*/
proc gammaii(a,p);
local i, j, g, l, m, n;
clear l;
g = zeros(rows(a),cols(a));
if rows(p) == rows(a) and cols(p) == cols(a);
l = 1;
m = rows(a);
n = cols(a);
elseif rows(p) == rows(a);
if cols(p) == 1;
l = 2;
m = rows(a);
n = cols(a);
elseif cols(a) == 1;
l = 3;
m = rows(p);
n = cols(p);
endif;
elseif cols(p) == cols(a);
if rows(p) == 1;
l = 4;
m = rows(a);
n = cols(a);
elseif rows(a) == 1;
l = 5;
m = rows(p);
n = cols(p);
endif;
endif;
if not l;
if not trapchk(1);
errorlog "GammaIncInv: arguments not conformable";
end;
else;
retp(error(0));
endif;
endif;
i = 1;
do until i > m;
j = 1;
do until j > n;
if l == 1;
g[i,j] = _ginvinc(a[i,j],p[i,j]);
elseif l == 2;
g[i,j] = _ginvinc(a[i,j],p[i]);
elseif l == 3;
g[i,j] = _ginvinc(a[i],p[i,j]);
elseif l == 4;
g[i,j] = _ginvinc(a[i,j],p[j]);
elseif l == 5;
g[i,j] = _ginvinc(a[j],p[i,j]);
endif;
j = j + 1;
endo;
i = i + 1;
endo;
retp(g);
endp;
proc _ginvinc(a,p);
local ln10;
local tol;
local euler;
local amax, b, d, g, h, q;
local s, t, u, r, w, y, z;
local c1, c2, c3, c4, c5, s2;
local qg, pn, qn, xn;
local am1, ap1, ap2, ap3, apn, rta;
local eps, sum;
local x,xmin,xmax;
if (p < 0 or p > 1) or a < 0;
if trapchk(1);
errorlog "GAMMII: arguments out of range";
end;
else;
retp(error(0));
endif;
endif;
ln10 = 2.3025850929940459;
tol = 1e-5;
euler = .577215664901533; /* Euler's Constant */
x = 0.0;
xmin = 2.2250738585072014e-308;
xmax = 1.7976931348623157e+308;
q = 1.0 - p;
if (p == 0.0);
retp(x);
endif;
if (q == 0.0);
retp(xmax);
endif;
if (a == 1.0);
if (q < 0.9);
retp(-ln(q));
else;
retp(-_gammaii_1(-p));
endif;
endif;
amax = 4e-11 / (__macheps * __macheps);
eps = 1e-10;
if (a <= 1.0);
g = gamma(a + 1.0);
qg = q * g;
if (qg == 0.0);
retp(xmax);
endif;
b = qg / a;
if (qg <= a * 0.6);
if (a < 0.3) and (b >= 0.35);
t = exp(-(b + euler));
u = t * exp(t);
xn = t * exp(u);
if (p > 0.5);
goto B;
else;
goto A;
endif;
endif;
if (b < 0.45);
if (b == 0.0);
retp(xmax);
endif;
y = -ln(b);
s = 0.5 - a + 0.5;
z = ln(y);
t = y - s * z;
if (b >= 0.15);
xn = y - s * ln(t) - ln(s / (t + 1.0) + 1.0);
goto B;
endif;
if (b > 0.01);
u = ((t + (3. - a) * 2.0) * t + (2.0 - a) * ( 3. - a))
/ ((t + (5. - a)) * t + 2.0);
xn = y - s * ln(t) - ln(u);
goto B;
endif;
c1 = -s * z;
c2 = -s * (c1 + 1.0);
c3 = s * ((c1 * 0.5 + (2.0 - a)) * c1 + (2.5 - a * 1.5));
c4 = -s * (((c1 / 3. + (2.5 - a * 1.5)) * c1 + ((a - 6.) *
a + 7.)) * c1 + ((a * 11. - 46) * a + 47.) / 6.);
c5 = -s * ((((-c1 / 4. + (a * 11. - 17.) / 6.) * c1 + ((a
* -3. + 13.) * a - 13.)) * c1 + (((a * 2.0 - 25.) * a +
72.) * a - 61.) * 0.5) * c1 + (((a * 25. - 195.) * a +
477.) * a - 379.) / 12.);
xn = (((c5 / y + c4) / y + c3) / y + c2) / y + c1 + y;
if (a > 1.0) or (b > 1e-28);
goto B;
endif;
retp(xn);
endif;
endif;
if (b * q > 1e-8);
if (p <= 0.9);
xn = exp(ln(p * g) / a);
else;
xn = exp((_gammaii_1(-q) + _gammaii_5(a)) / a);
endif;
else;
xn = exp(-(q / a + euler));
endif;
t = 0.5 - xn / (a + 1.0) + 0.5;
xn = xn / t;
if (p > 0.5);
goto B;
else;
goto A;
endif;
else;
b = 0;
endif;
if (q <= 0.5);
w = ln(q);
else;
w = ln(p);
endif;
t = sqrt(w * -2.);
s = t - (((.213623493715853 * t + 4.28342155967104) * t +
11.6616720288968) * t + 3.31125922108741) / ((((0.036117081018842
* t + 1.27364489782223) * t + 6.40691597760039) * t +
6.61053765625462) * t + 1.0);
if (q > 0.5);
s = -s;
endif;
rta = sqrt(a);
s2 = s * s;
xn = a + s * rta + (s2 - 1.0) / 3. + s * (s2 - 7.) / (rta * 36.) -
((s2 * 3. + 7.) * s2 - 16.) / (a * 810.) + s * ((s2 * 9. + 256.) *
s2 - 433.) / (a * 38880. * rta);
xn = maxc(xn | 0);
if (a >= 500);
x = xn;
if (abs(1.0 - x / a) <= 1e-6);
retp(x);
endif;
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?