Doxygen Source Code Documentation
cdf_23.c File Reference
#include "cdflib.h"
Go to the source code of this file.
Defines | |
#define | tol (1.0e-8) |
#define | atol (1.0e-50) |
#define | inf 1.0e300 |
Functions | |
void | cdfpoi (int *which, double *p, double *q, double *s, double *xlam, int *status, double *bound) |
Define Documentation
|
|
|
|
|
|
Function Documentation
|
Definition at line 2 of file cdf_23.c. References cumpoi(), dinvr(), dstinv(), p, q, and spmpar(). Referenced by poisson_p2t(), and poisson_t2p().
00025 : 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 */ |