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