|
Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
cdf_18.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 | zero (1.0e-300) |
#define | inf 1.0e300 |
Functions |
void | cdff (int *which, double *p, double *q, double *f, double *dfn, double *dfd, int *status, double *bound) |
Define Documentation
Function Documentation
void cdff |
( |
int * |
which, |
|
|
double * |
p, |
|
|
double * |
q, |
|
|
double * |
f, |
|
|
double * |
dfn, |
|
|
double * |
dfd, |
|
|
int * |
status, |
|
|
double * |
bound |
|
) |
|
|
|
Definition at line 2 of file cdf_18.c.
References cumf(), dinvr(), dstinv(), p, q, and spmpar().
Referenced by fstat_p2t(), fstat_t2p(), fstat_t2pp(), and identify_repeats().
00025 : 1..4
00026 iwhich = 1 : Calculate P and Q from F,DFN and DFD
00027 iwhich = 2 : Calculate F from P,Q,DFN and DFD
00028 iwhich = 3 : Calculate DFN from P,Q,F and DFD
00029 iwhich = 4 : Calculate DFD from P,Q,F and DFN
00030
00031 P <--> The integral from 0 to F of the f-density.
00032 Input range: [0,1].
00033
00034 Q <--> 1-P.
00035 Input range: (0, 1].
00036 P + Q = 1.0.
00037
00038 F <--> Upper limit of integration of the f-density.
00039 Input range: [0, +infinity).
00040 Search range: [0,1E300]
00041
00042 DFN < --> Degrees of freedom of the numerator sum of squares.
00043 Input range: (0, +infinity).
00044 Search range: [ 1E-300, 1E300]
00045
00046 DFD < --> Degrees of freedom of the denominator sum of squares.
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
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.6.2 of Abramowitz and Stegun, Handbook of
00072 Mathematical Functions (1966) is used to reduce the computation
00073 of the cumulative distribution function for the F variate to
00074 that of an incomplete beta.
00075
00076 Computation of other parameters involve a seach for a value that
00077 produces the desired value of P. The search relies on the
00078 monotinicity of P with the other parameter.
00079
00080 WARNING
00081
00082 The value of the cumulative F distribution is not necessarily
00083 monotone in either degrees of freedom. There thus may be two
00084 values that provide a given CDF value. This routine assumes
00085 monotonicity and will find an arbitrary one of the two values.
00086
00087 **********************************************************************/
00088 {
00089 #define tol (1.0e-8)
00090 #define atol (1.0e-50)
00091 #define zero (1.0e-300)
00092 #define inf 1.0e300
00093 static int K1 = 1;
00094 static double K2 = 0.0e0;
00095 static double K4 = 0.5e0;
00096 static double K5 = 5.0e0;
00097 static double pq,fx,cum,ccum;
00098 static unsigned long qhi,qleft,qporq;
00099 static double T3,T6,T7,T8,T9,T10,T11,T12,T13,T14,T15;
00100
00101
00102
00103
00104
00105
00106
00107 if(!(*which < 1 || *which > 4)) goto S30;
00108 if(!(*which < 1)) goto S10;
00109 *bound = 1.0e0;
00110 goto S20;
00111 S10:
00112 *bound = 4.0e0;
00113 S20:
00114 *status = -1;
00115 return;
00116 S30:
00117 if(*which == 1) goto S70;
00118
00119
00120
00121 if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
00122 if(!(*p < 0.0e0)) goto S40;
00123 *bound = 0.0e0;
00124 goto S50;
00125 S40:
00126 *bound = 1.0e0;
00127 S50:
00128 *status = -2;
00129 return;
00130 S70:
00131 S60:
00132 if(*which == 1) goto S110;
00133
00134
00135
00136 if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
00137 if(!(*q <= 0.0e0)) goto S80;
00138 *bound = 0.0e0;
00139 goto S90;
00140 S80:
00141 *bound = 1.0e0;
00142 S90:
00143 *status = -3;
00144 return;
00145 S110:
00146 S100:
00147 if(*which == 2) goto S130;
00148
00149
00150
00151 if(!(*f < 0.0e0)) goto S120;
00152 *bound = 0.0e0;
00153 *status = -4;
00154 return;
00155 S130:
00156 S120:
00157 if(*which == 3) goto S150;
00158
00159
00160
00161 if(!(*dfn <= 0.0e0)) goto S140;
00162 *bound = 0.0e0;
00163 *status = -5;
00164 return;
00165 S150:
00166 S140:
00167 if(*which == 4) goto S170;
00168
00169
00170
00171 if(!(*dfd <= 0.0e0)) goto S160;
00172 *bound = 0.0e0;
00173 *status = -6;
00174 return;
00175 S170:
00176 S160:
00177 if(*which == 1) goto S210;
00178
00179
00180
00181 pq = *p+*q;
00182 if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S200;
00183 if(!(pq < 0.0e0)) goto S180;
00184 *bound = 0.0e0;
00185 goto S190;
00186 S180:
00187 *bound = 1.0e0;
00188 S190:
00189 *status = 3;
00190 return;
00191 S210:
00192 S200:
00193 if(!(*which == 1)) qporq = *p <= *q;
00194
00195
00196
00197
00198 if(1 == *which) {
00199
00200
00201
00202 cumf(f,dfn,dfd,p,q);
00203 *status = 0;
00204 }
00205 else if(2 == *which) {
00206
00207
00208
00209 *f = 5.0e0;
00210 T3 = inf;
00211 T6 = atol;
00212 T7 = tol;
00213 dstinv(&K2,&T3,&K4,&K4,&K5,&T6,&T7);
00214 *status = 0;
00215 dinvr(status,f,&fx,&qleft,&qhi);
00216 S220:
00217 if(!(*status == 1)) goto S250;
00218 cumf(f,dfn,dfd,&cum,&ccum);
00219 if(!qporq) goto S230;
00220 fx = cum-*p;
00221 goto S240;
00222 S230:
00223 fx = ccum-*q;
00224 S240:
00225 dinvr(status,f,&fx,&qleft,&qhi);
00226 goto S220;
00227 S250:
00228 if(!(*status == -1)) goto S280;
00229 if(!qleft) goto S260;
00230 *status = 1;
00231 *bound = 0.0e0;
00232 goto S270;
00233 S260:
00234 *status = 2;
00235 *bound = inf;
00236 S280:
00237 S270:
00238 ;
00239 }
00240 else if(3 == *which) {
00241
00242
00243
00244 *dfn = 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,dfn,&fx,&qleft,&qhi);
00252 S290:
00253 if(!(*status == 1)) goto S320;
00254 cumf(f,dfn,dfd,&cum,&ccum);
00255 if(!qporq) goto S300;
00256 fx = cum-*p;
00257 goto S310;
00258 S300:
00259 fx = ccum-*q;
00260 S310:
00261 dinvr(status,dfn,&fx,&qleft,&qhi);
00262 goto S290;
00263 S320:
00264 if(!(*status == -1)) goto S350;
00265 if(!qleft) goto S330;
00266 *status = 1;
00267 *bound = zero;
00268 goto S340;
00269 S330:
00270 *status = 2;
00271 *bound = inf;
00272 S350:
00273 S340:
00274 ;
00275 }
00276 else if(4 == *which) {
00277
00278
00279
00280 *dfd = 5.0e0;
00281 T12 = zero;
00282 T13 = inf;
00283 T14 = atol;
00284 T15 = tol;
00285 dstinv(&T12,&T13,&K4,&K4,&K5,&T14,&T15);
00286 *status = 0;
00287 dinvr(status,dfd,&fx,&qleft,&qhi);
00288 S360:
00289 if(!(*status == 1)) goto S390;
00290 cumf(f,dfn,dfd,&cum,&ccum);
00291 if(!qporq) goto S370;
00292 fx = cum-*p;
00293 goto S380;
00294 S370:
00295 fx = ccum-*q;
00296 S380:
00297 dinvr(status,dfd,&fx,&qleft,&qhi);
00298 goto S360;
00299 S390:
00300 if(!(*status == -1)) goto S420;
00301 if(!qleft) goto S400;
00302 *status = 1;
00303 *bound = zero;
00304 goto S410;
00305 S400:
00306 *status = 2;
00307 *bound = inf;
00308 S410:
00309 ;
00310 }
00311 S420:
00312 return;
00313 #undef tol
00314 #undef atol
00315 #undef zero
00316 #undef inf
00317 }
|
|