Math-CDF
view release on metacpan or search on metacpan
cdflib/dcdflib.c view on Meta::CPAN
d = *b+(*a-0.5e0);
S20:
/*
SET SN = (1 - X**N)/(1 - X)
*/
x2 = x*x;
s3 = 1.0e0+(x+x2);
s5 = 1.0e0+(x+x2*s3);
s7 = 1.0e0+(x+x2*s5);
s9 = 1.0e0+(x+x2*s7);
s11 = 1.0e0+(x+x2*s9);
/*
SET W = DEL(B) - DEL(A + B)
*/
t = pow(1.0e0/ *b,2.0);
w = ((((c5*s11*t+c4*s9)*t+c3*s7)*t+c2*s5)*t+c1*s3)*t+c0;
w *= (c/ *b);
/*
COMBINE THE RESULTS
*/
T1 = *a/ *b;
u = d*alnrel(&T1);
v = *a*(log(*b)-1.0e0);
if(u <= v) goto S30;
algdiv = w-v-u;
return algdiv;
S30:
algdiv = w-u-v;
return algdiv;
}
double alngam(double *x)
/*
**********************************************************************
double alngam(double *x)
double precision LN of the GAMma function
Function
Returns the natural logarithm of GAMMA(X).
Arguments
X --> value at which scaled log gamma is to be returned
X is DOUBLE PRECISION
Method
If X .le. 6.0, then use recursion to get X below 3
then apply rational approximation number 5236 of
Hart et al, Computer Approximations, John Wiley and
Sons, NY, 1968.
If X .gt. 6.0, then use recursion to get X to at least 12 and
then use formula 5423 of the same source.
**********************************************************************
*/
{
#define hln2pi 0.91893853320467274178e0
static double coef[5] = {
0.83333333333333023564e-1,-0.27777777768818808e-2,0.79365006754279e-3,
-0.594997310889e-3,0.8065880899e-3
};
static double scoefd[4] = {
0.62003838007126989331e2,0.9822521104713994894e1,-0.8906016659497461257e1,
0.1000000000000000000e1
};
static double scoefn[9] = {
0.62003838007127258804e2,0.36036772530024836321e2,0.20782472531792126786e2,
0.6338067999387272343e1,0.215994312846059073e1,0.3980671310203570498e0,
0.1093115956710439502e0,0.92381945590275995e-2,0.29737866448101651e-2
};
static int K1 = 9;
static int K3 = 4;
static int K5 = 5;
static double alngam,offset,prod,xx;
static int i,n;
static double T2,T4,T6;
/*
..
.. Executable Statements ..
*/
if(!(*x <= 6.0e0)) goto S70;
prod = 1.0e0;
xx = *x;
if(!(*x > 3.0e0)) goto S30;
S10:
if(!(xx > 3.0e0)) goto S20;
xx -= 1.0e0;
prod *= xx;
goto S10;
S30:
S20:
if(!(*x < 2.0e0)) goto S60;
S40:
if(!(xx < 2.0e0)) goto S50;
prod /= xx;
xx += 1.0e0;
goto S40;
S60:
S50:
T2 = xx-2.0e0;
T4 = xx-2.0e0;
alngam = devlpl(scoefn,&K1,&T2)/devlpl(scoefd,&K3,&T4);
/*
COMPUTE RATIONAL APPROXIMATION TO GAMMA(X)
*/
alngam *= prod;
alngam = log(alngam);
goto S110;
S70:
offset = hln2pi;
/*
IF NECESSARY MAKE X AT LEAST 12 AND CARRY CORRECTION IN OFFSET
cdflib/dcdflib.c view on Meta::CPAN
iwhich = 1 : Calculate P and Q from X,MEAN and SD
iwhich = 2 : Calculate X from P,Q,MEAN and SD
iwhich = 3 : Calculate MEAN from P,Q,X and SD
iwhich = 4 : Calculate SD from P,Q,X and MEAN
P <--> The integral from -infinity to X of the normal density.
Input range: (0,1].
Q <--> 1-P.
Input range: (0, 1].
P + Q = 1.0.
X < --> Upper limit of integration of the normal-density.
Input range: ( -infinity, +infinity)
MEAN <--> The mean of the normal density.
Input range: (-infinity, +infinity)
SD <--> Standard Deviation of the normal density.
Input range: (0, +infinity).
STATUS <-- 0 if calculation completed correctly
-I if input parameter number I is out of range
1 if answer appears to be lower than lowest
search bound
2 if answer appears to be higher than greatest
search bound
3 if P + Q .ne. 1
BOUND <-- Undefined if STATUS is 0
Bound exceeded by parameter number I if STATUS
is negative.
Lower search bound if STATUS is 1.
Upper search bound if STATUS is 2.
Method
A slightly modified version of ANORM from
Cody, W.D. (1993). "ALGORITHM 715: SPECFUN - A Portabel FORTRAN
Package of Special Function Routines and Test Drivers"
acm Transactions on Mathematical Software. 19, 22-32.
is used to calulate the cumulative standard normal distribution.
The rational functions from pages 90-95 of Kennedy and Gentle,
Statistical Computing, Marcel Dekker, NY, 1980 are used as
starting values to Newton's Iterations which compute the inverse
standard normal. Therefore no searches are necessary for any
parameter.
For X < -15, the asymptotic expansion for the normal is used as
the starting value in finding the inverse standard normal.
This is formula 26.2.12 of Abramowitz and Stegun.
Note
The normal density is proportional to
exp( - 0.5 * (( X - MEAN)/SD)**2)
**********************************************************************/
{
static int K1 = 1;
static double z,pq;
/*
..
.. Executable Statements ..
*/
/*
Check arguments
*/
*status = 0;
if(!(*which < 1 || *which > 4)) goto S30;
if(!(*which < 1)) goto S10;
*bound = 1.0e0;
goto S20;
S10:
*bound = 4.0e0;
S20:
*status = -1;
return;
S30:
if(*which == 1) goto S70;
/*
P
*/
if(!(*p <= 0.0e0 || *p > 1.0e0)) goto S60;
if(!(*p <= 0.0e0)) goto S40;
*bound = 0.0e0;
goto S50;
S40:
*bound = 1.0e0;
S50:
*status = -2;
return;
S70:
S60:
if(*which == 1) goto S110;
/*
Q
*/
if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
if(!(*q <= 0.0e0)) goto S80;
*bound = 0.0e0;
goto S90;
S80:
*bound = 1.0e0;
S90:
*status = -3;
return;
S110:
S100:
cdflib/dcdflib.c view on Meta::CPAN
void cdftnc(int *which,double *p,double *q,double *t,double *df,
double *pnonc,int *status,double *bound)
Cumulative Distribution Function
Non-Central T distribution
Function
Calculates any one parameter of the noncentral t distribution give
values for the others.
Arguments
WHICH --> Integer indicating which argument
values is to be calculated from the others.
Legal range: 1..3
iwhich = 1 : Calculate P and Q from T,DF,PNONC
iwhich = 2 : Calculate T from P,Q,DF,PNONC
iwhich = 3 : Calculate DF from P,Q,T
iwhich = 4 : Calculate PNONC from P,Q,DF,T
P <--> The integral from -infinity to t of the noncentral t-den
Input range: (0,1].
Q <--> 1-P.
Input range: (0, 1].
P + Q = 1.0.
T <--> Upper limit of integration of the noncentral t-density.
Input range: ( -infinity, +infinity).
Search range: [ -1E100, 1E100 ]
DF <--> Degrees of freedom of the noncentral t-distribution.
Input range: (0 , +infinity).
Search range: [1e-100, 1E10]
PNONC <--> Noncentrality parameter of the noncentral t-distributio
Input range: [-infinity , +infinity).
Search range: [-1e4, 1E4]
STATUS <-- 0 if calculation completed correctly
-I if input parameter number I is out of range
1 if answer appears to be lower than lowest
search bound
2 if answer appears to be higher than greatest
search bound
3 if P + Q .ne. 1
BOUND <-- Undefined if STATUS is 0
Bound exceeded by parameter number I if STATUS
is negative.
Lower search bound if STATUS is 1.
Upper search bound if STATUS is 2.
Method
Upper tail of the cumulative noncentral t is calculated usin
formulae from page 532 of Johnson, Kotz, Balakrishnan, Coninuou
Univariate Distributions, Vol 2, 2nd Edition. Wiley (1995)
Computation of other parameters involve a seach for a value that
produces the desired value of P. The search relies on the
monotinicity of P with the other parameter.
**********************************************************************/
{
#define tent4 1.0e4
#define tol 1.0e-8
#define atol 1.0e-50
#define zero 1.0e-100
#define one ( 1.0e0 - 1.0e-16 )
#define inf 1.0e100
static double K3 = 0.5e0;
static double K4 = 5.0e0;
static double ccum,cum,fx;
static unsigned long qhi,qleft;
static double T1,T2,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14;
/*
..
.. Executable Statements ..
*/
if(!(*which < 1 || *which > 4)) goto S30;
if(!(*which < 1)) goto S10;
*bound = 1.0e0;
goto S20;
S10:
*bound = 5.0e0;
S20:
*status = -1;
return;
S30:
if(*which == 1) goto S70;
if(!(*p < 0.0e0 || *p > one)) goto S60;
if(!(*p < 0.0e0)) goto S40;
*bound = 0.0e0;
goto S50;
S40:
*bound = one;
S50:
*status = -2;
return;
S70:
S60:
if(*which == 3) goto S90;
if(!(*df <= 0.0e0)) goto S80;
*bound = 0.0e0;
*status = -5;
return;
S90:
S80:
if(*which == 4) goto S100;
S100:
if(1 == *which) {
cumtnc(t,df,pnonc,p,q);
*status = 0;
}
else if(2 == *which) {
*t = 5.0e0;
cdflib/dcdflib.c view on Meta::CPAN
CUM <-- Cumulative chi-square distribution.
CUM is DOUBLE PRECISIO
CCUM <-- Compliment of Cumulative chi-square distribution.
CCUM is DOUBLE PRECISI
Method
Calls incomplete gamma function (CUMGAM)
**********************************************************************
*/
{
static double a,xx;
/*
..
.. Executable Statements ..
*/
a = *df*0.5e0;
xx = *x*0.5e0;
cumgam(&xx,&a,cum,ccum);
return;
}
void cumchn(double *x,double *df,double *pnonc,double *cum,
double *ccum)
/**********************************************************************
void cumchn(double *x,double *df,double *pnonc,double *cum,
double *ccum)
CUMulative of the Non-central CHi-square distribution
Function
Calculates the cumulative non-central chi-square
distribution, i.e., the probability that a random variable
which follows the non-central chi-square distribution, with
non-centrality parameter PNONC and continuous degrees of
freedom DF, is less than or equal to X.
Arguments
X --> Upper limit of integration of the non-central
chi-square distribution.
DF --> Degrees of freedom of the non-central
chi-square distribution.
PNONC --> Non-centrality parameter of the non-central
chi-square distribution.
CUM <-- Cumulative non-central chi-square distribution.
CCUM <-- Compliment of Cumulative non-central chi-square distribut
Method
Uses formula 26.4.25 of Abramowitz and Stegun, Handbook of
Mathematical Functions, US NBS (1966) to calculate the
non-central chi-square.
Variables
EPS --- Convergence criterion. The sum stops when a
term is less than EPS*SUM.
CCUM <-- Compliment of Cumulative non-central
chi-square distribution.
**********************************************************************/
{
#define dg(i) (*df + 2.0e0 * (double)(i))
#define qsmall(xx) (int)(sum < 1.0e-20 || (xx) < eps * sum)
static double eps = 1.0e-5;
static double adj,centaj,centwt,chid2,dfd2,lcntaj,lcntwt,lfact,pcent,pterm,sum,
sumadj,term,wt,xnonc;
static int i,icent;
static double T1,T2,T3;
/*
..
.. Executable Statements ..
*/
if(!(*x <= 0.0e0)) goto S10;
*cum = 0.0e0;
*ccum = 1.0e0;
return;
S10:
if(!(*pnonc <= 1.0e-10 )) goto S20;
/*
When non-centrality parameter is (essentially) zero,
use cumulative chi-square distribution
*/
cumchi(x,df,cum,ccum);
return;
S20:
xnonc = *pnonc / 2.0e0;
/*
***********************************************************************
The following code calcualtes the weight, chi-square, and
adjustment term for the central term in the infinite series.
The central term is the one in which the poisson weight is
greatest. The adjustment term is the amount that must
be subtracted from the chi-square to move up two degrees
of freedom.
***********************************************************************
*/
icent = fifidint(xnonc);
if(icent == 0) icent = 1;
chid2 = *x / 2.0e0;
/*
Calculate central weight term
*/
T1 = (double)(icent + 1);
lfact = alngam(&T1);
lcntwt = -xnonc + (double)icent * log(xnonc) - lfact;
centwt = exp(lcntwt);
/*
Calculate central chi-square
cdflib/dcdflib.c view on Meta::CPAN
xden = xsq;
for(i=0; i<4; i++) {
xnum = (xnum+p[i])*xsq;
xden = (xden+q[i])*xsq;
}
*result = xsq*(xnum+p[4])/(xden+q[4]);
*result = (sqrpi-*result)/y;
xsq = fifdint(x*sixten)/sixten;
del = (x-xsq)*(x+xsq);
*result = exp(-(xsq*xsq*half))*exp(-(del*half))**result;
*ccum = one-*result;
if(x > zero) {
temp = *result;
*result = *ccum;
*ccum = temp;
}
}
if(*result < min) *result = 0.0e0;
/*
------------------------------------------------------------------
Fix up for negative argument, erf, etc.
------------------------------------------------------------------
----------Last card of ANORM ----------
*/
if(*ccum < min) *ccum = 0.0e0;
}
void cumpoi(double *s,double *xlam,double *cum,double *ccum)
/*
**********************************************************************
void cumpoi(double *s,double *xlam,double *cum,double *ccum)
CUMulative POIsson distribution
Function
Returns the probability of S or fewer events in a Poisson
distribution with mean XLAM.
Arguments
S --> Upper limit of cumulation of the Poisson.
S is DOUBLE PRECISION
XLAM --> Mean of the Poisson distribution.
XLAM is DOUBLE PRECIS
CUM <-- Cumulative poisson distribution.
CUM is DOUBLE PRECISION
CCUM <-- Compliment of Cumulative poisson distribution.
CCUM is DOUBLE PRECIS
Method
Uses formula 26.4.21 of Abramowitz and Stegun, Handbook of
Mathematical Functions to reduce the cumulative Poisson to
the cumulative chi-square distribution.
**********************************************************************
*/
{
static double chi,df;
/*
..
.. Executable Statements ..
*/
df = 2.0e0*(*s+1.0e0);
chi = 2.0e0**xlam;
cumchi(&chi,&df,ccum,cum);
return;
}
void cumt(double *t,double *df,double *cum,double *ccum)
/*
**********************************************************************
void cumt(double *t,double *df,double *cum,double *ccum)
CUMulative T-distribution
Function
Computes the integral from -infinity to T of the t-density.
Arguments
T --> Upper limit of integration of the t-density.
T is DOUBLE PRECISION
DF --> Degrees of freedom of the t-distribution.
DF is DOUBLE PRECISIO
CUM <-- Cumulative t-distribution.
CCUM is DOUBLE PRECIS
CCUM <-- Compliment of Cumulative t-distribution.
CCUM is DOUBLE PRECIS
Method
Formula 26.5.27 of Abramowitz and Stegun, Handbook of
Mathematical Functions is used to reduce the t-distribution
to an incomplete beta.
**********************************************************************
*/
{
static double K2 = 0.5e0;
static double xx,a,oma,tt,yy,dfptt,T1;
/*
..
.. Executable Statements ..
*/
tt = *t**t;
dfptt = *df+tt;
xx = *df/dfptt;
yy = tt/dfptt;
T1 = 0.5e0**df;
cumbet(&xx,&yy,&T1,&K2,&a,&oma);
if(!(*t <= 0.0e0)) goto S10;
*cum = 0.5e0*a;
*ccum = oma+*cum;
goto S20;
S10:
*ccum = 0.5e0*a;
*cum = oma+*ccum;
S20:
return;
}
void cumtnc(double *t,double *df,double *pnonc,double *cum,
double *ccum)
/**********************************************************************
void cumtnc(double *t,double *df,double *pnonc,double *cum,
double *ccum)
CUMulative Non-Central T-distribution
Function
Computes the integral from -infinity to T of the non-central
t-density.
Arguments
T --> Upper limit of integration of the non-central t-density.
DF --> Degrees of freedom of the non-central t-distribution.
PNONC --> Non-centrality parameter of the non-central t distibutio
CUM <-- Cumulative t-distribution.
CCUM <-- Compliment of Cumulative t-distribution.
Method
Upper tail of the cumulative noncentral t using
formulae from page 532 of Johnson, Kotz, Balakrishnan, Coninuous
Univariate Distributions, Vol 2, 2nd Edition. Wiley (1995)
This implementation starts the calculation at i = lambda,
which is near the largest Di. It then sums forward and backward.
**********************************************************************/
{
#define one 1.0e0
#define zero 0.0e0
#define half 0.5e0
#define two 2.0e0
#define onep5 1.5e0
#define conv 1.0e-7
#define tiny 1.0e-10
static double alghdf,b,bb,bbcent,bcent,cent,d,dcent,dpnonc,dum1,dum2,e,ecent,
halfdf,lambda,lnomx,lnx,omx,pnonc2,s,scent,ss,sscent,t2,term,tt,twoi,x,xi,
xlnd,xlne;
static int ierr;
static unsigned long qrevs;
static double T1,T2,T3,T4,T5,T6,T7,T8,T9,T10;
/*
..
.. Executable Statements ..
*/
/*
Case pnonc essentially zero
*/
if(fabs(*pnonc) <= tiny) {
cumt(t,df,cum,ccum);
return;
}
qrevs = *t < zero;
if(qrevs) {
tt = -*t;
dpnonc = -*pnonc;
}
else {
tt = *t;
dpnonc = *pnonc;
}
pnonc2 = dpnonc * dpnonc;
t2 = tt * tt;
if(fabs(tt) <= tiny) {
T1 = -*pnonc;
cumnor(&T1,cum,ccum);
return;
}
lambda = half * pnonc2;
x = *df / (*df + t2);
omx = one - x;
lnx = log(x);
lnomx = log(omx);
halfdf = half * *df;
alghdf = gamln(&halfdf);
/*
******************** Case i = lambda
*/
cent = fifidint(lambda);
if(cent < one) cent = one;
/*
Compute d=T(2i) in log space and offset by exp(-lambda)
( run in 1.270 second using v1.01-cache-2.11-cpan-39bf76dae61 )