00001 #include "cdflib.h"
00002 void cdfchi(int *which,double *p,double *q,double *x,double *df,
00003             int *status,double *bound)
00004 /**********************************************************************
00005 
00006       void cdfchi(int *which,double *p,double *q,double *x,double *df,
00007             int *status,double *bound)
00008 
00009                Cumulative Distribution Function
00010                CHI-Square distribution
00011 
00012 
00013                               Function
00014 
00015 
00016      Calculates any one parameter of the chi-square
00017      distribution given values for the others.
00018 
00019 
00020                               Arguments
00021 
00022 
00023      WHICH --> Integer indicating which of the next three argument
00024                values is to be calculated from the others.
00025                Legal range: 1..3
00026                iwhich = 1 : Calculate P and Q from X and DF
00027                iwhich = 2 : Calculate X from P,Q and DF
00028                iwhich = 3 : Calculate DF from P,Q and X
00029 
00030      P <--> The integral from 0 to X of the chi-square
00031             distribution.
00032             Input range: [0, 1].
00033 
00034      Q <--> 1-P.
00035             Input range: (0, 1].
00036             P + Q = 1.0.
00037 
00038      X <--> Upper limit of integration of the non-central
00039             chi-square distribution.
00040             Input range: [0, +infinity).
00041             Search range: [0,1E300]
00042 
00043      DF <--> Degrees of freedom of the
00044              chi-square distribution.
00045              Input range: (0, +infinity).
00046              Search range: [ 1E-300, 1E300]
00047 
00048      STATUS <-- 0 if calculation completed correctly
00049                -I if input parameter number I is out of range
00050                 1 if answer appears to be lower than lowest
00051                   search bound
00052                 2 if answer appears to be higher than greatest
00053                   search bound
00054                 3 if P + Q .ne. 1
00055                10 indicates error returned from cumgam.  See
00056                   references in cdfgam
00057 
00058      BOUND <-- Undefined if STATUS is 0
00059 
00060                Bound exceeded by parameter number I if STATUS
00061                is negative.
00062 
00063                Lower search bound if STATUS is 1.
00064 
00065                Upper search bound if STATUS is 2.
00066 
00067 
00068                               Method
00069 
00070 
00071      Formula    26.4.19   of Abramowitz  and     Stegun, Handbook  of
00072      Mathematical Functions   (1966) is used   to reduce the chisqure
00073      distribution to the incomplete distribution.
00074 
00075      Computation of other parameters involve a seach for a value that
00076      produces  the desired  value  of P.   The search relies  on  the
00077      monotinicity of P with the other parameter.
00078 
00079 **********************************************************************/
00080 {
00081 #define tol (1.0e-8)
00082 #define atol (1.0e-50)
00083 #define zero (1.0e-300)
00084 #define inf 1.0e300
00085 static int K1 = 1;
00086 static double K2 = 0.0e0;
00087 static double K4 = 0.5e0;
00088 static double K5 = 5.0e0;
00089 static double fx,cum,ccum,pq,porq;
00090 static unsigned long qhi,qleft,qporq;
00091 static double T3,T6,T7,T8,T9,T10,T11;
00092 /*
00093      ..
00094      .. Executable Statements ..
00095 */
00096 /*
00097      Check arguments
00098 */
00099     if(!(*which < 1 || *which > 3)) goto S30;
00100     if(!(*which < 1)) goto S10;
00101     *bound = 1.0e0;
00102     goto S20;
00103 S10:
00104     *bound = 3.0e0;
00105 S20:
00106     *status = -1;
00107     return;
00108 S30:
00109     if(*which == 1) goto S70;
00110 /*
00111      P
00112 */
00113     if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
00114     if(!(*p < 0.0e0)) goto S40;
00115     *bound = 0.0e0;
00116     goto S50;
00117 S40:
00118     *bound = 1.0e0;
00119 S50:
00120     *status = -2;
00121     return;
00122 S70:
00123 S60:
00124     if(*which == 1) goto S110;
00125 /*
00126      Q
00127 */
00128     if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
00129     if(!(*q <= 0.0e0)) goto S80;
00130     *bound = 0.0e0;
00131     goto S90;
00132 S80:
00133     *bound = 1.0e0;
00134 S90:
00135     *status = -3;
00136     return;
00137 S110:
00138 S100:
00139     if(*which == 2) goto S130;
00140 /*
00141      X
00142 */
00143     if(!(*x < 0.0e0)) goto S120;
00144     *bound = 0.0e0;
00145     *status = -4;
00146     return;
00147 S130:
00148 S120:
00149     if(*which == 3) goto S150;
00150 /*
00151      DF
00152 */
00153     if(!(*df <= 0.0e0)) goto S140;
00154     *bound = 0.0e0;
00155     *status = -5;
00156     return;
00157 S150:
00158 S140:
00159     if(*which == 1) goto S190;
00160 /*
00161      P + Q
00162 */
00163     pq = *p+*q;
00164     if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S180;
00165     if(!(pq < 0.0e0)) goto S160;
00166     *bound = 0.0e0;
00167     goto S170;
00168 S160:
00169     *bound = 1.0e0;
00170 S170:
00171     *status = 3;
00172     return;
00173 S190:
00174 S180:
00175     if(*which == 1) goto S220;
00176 /*
00177      Select the minimum of P or Q
00178 */
00179     qporq = *p <= *q;
00180     if(!qporq) goto S200;
00181     porq = *p;
00182     goto S210;
00183 S200:
00184     porq = *q;
00185 S220:
00186 S210:
00187 /*
00188      Calculate ANSWERS
00189 */
00190     if(1 == *which) {
00191 /*
00192      Calculating P and Q
00193 */
00194         *status = 0;
00195         cumchi(x,df,p,q);
00196         if(porq > 1.5e0) {
00197             *status = 10;
00198             return;
00199         }
00200     }
00201     else if(2 == *which) {
00202 /*
00203      Calculating X
00204 */
00205         *x = 5.0e0;
00206         T3 = inf;
00207         T6 = atol;
00208         T7 = tol;
00209         dstinv(&K2,&T3,&K4,&K4,&K5,&T6,&T7);
00210         *status = 0;
00211         dinvr(status,x,&fx,&qleft,&qhi);
00212 S230:
00213         if(!(*status == 1)) goto S270;
00214         cumchi(x,df,&cum,&ccum);
00215         if(!qporq) goto S240;
00216         fx = cum-*p;
00217         goto S250;
00218 S240:
00219         fx = ccum-*q;
00220 S250:
00221         if(!(fx+porq > 1.5e0)) goto S260;
00222         *status = 10;
00223         return;
00224 S260:
00225         dinvr(status,x,&fx,&qleft,&qhi);
00226         goto S230;
00227 S270:
00228         if(!(*status == -1)) goto S300;
00229         if(!qleft) goto S280;
00230         *status = 1;
00231         *bound = 0.0e0;
00232         goto S290;
00233 S280:
00234         *status = 2;
00235         *bound = inf;
00236 S300:
00237 S290:
00238         ;
00239     }
00240     else if(3 == *which) {
00241 /*
00242      Calculating DF
00243 */
00244         *df = 5.0e0;
00245         T8 = zero;
00246         T9 = inf;
00247         T10 = atol;
00248         T11 = tol;
00249         dstinv(&T8,&T9,&K4,&K4,&K5,&T10,&T11);
00250         *status = 0;
00251         dinvr(status,df,&fx,&qleft,&qhi);
00252 S310:
00253         if(!(*status == 1)) goto S350;
00254         cumchi(x,df,&cum,&ccum);
00255         if(!qporq) goto S320;
00256         fx = cum-*p;
00257         goto S330;
00258 S320:
00259         fx = ccum-*q;
00260 S330:
00261         if(!(fx+porq > 1.5e0)) goto S340;
00262         *status = 10;
00263         return;
00264 S340:
00265         dinvr(status,df,&fx,&qleft,&qhi);
00266         goto S310;
00267 S350:
00268         if(!(*status == -1)) goto S380;
00269         if(!qleft) goto S360;
00270         *status = 1;
00271         *bound = zero;
00272         goto S370;
00273 S360:
00274         *status = 2;
00275         *bound = inf;
00276 S370:
00277         ;
00278     }
00279 S380:
00280     return;
00281 #undef tol
00282 #undef atol
00283 #undef zero
00284 #undef inf
00285 } /* END */