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