|
|
Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
cdf_14.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 |
| #define | one 1.0e0 |
Functions |
| void | cdfbet (int *which, double *p, double *q, double *x, double *y, double *a, double *b, int *status, double *bound) |
Define Documentation
Function Documentation
| void cdfbet |
( |
int * |
which, |
|
|
double * |
p, |
|
|
double * |
q, |
|
|
double * |
x, |
|
|
double * |
y, |
|
|
double * |
a, |
|
|
double * |
b, |
|
|
int * |
status, |
|
|
double * |
bound |
|
) |
|
|
|
|
Definition at line 2 of file cdf_14.c.
References a, cumbet(), dinvr(), dstinv(), dstzr(), dzror(), p, q, and spmpar().
Referenced by beta_p2t(), and beta_t2p().
00025 : 1..4
00026 iwhich = 1 : Calculate P and Q from X,Y,A and B
00027 iwhich = 2 : Calculate X and Y from P,Q,A and B
00028 iwhich = 3 : Calculate A from P,Q,X,Y and B
00029 iwhich = 4 : Calculate B from P,Q,X,Y and A
00030
00031 P <--> The integral from 0 to X of the chi-square
00032 distribution.
00033 Input range: [0, 1].
00034
00035 Q <--> 1-P.
00036 Input range: [0, 1].
00037 P + Q = 1.0.
00038
00039 X <--> Upper limit of integration of beta density.
00040 Input range: [0,1].
00041 Search range: [0,1]
00042
00043 Y <--> 1-X.
00044 Input range: [0,1].
00045 Search range: [0,1]
00046 X + Y = 1.0.
00047
00048 A <--> The first parameter of the beta density.
00049 Input range: (0, +infinity).
00050 Search range: [1D-300,1D300]
00051
00052 B <--> The second parameter of the beta density.
00053 Input range: (0, +infinity).
00054 Search range: [1D-300,1D300]
00055
00056 STATUS <-- 0 if calculation completed correctly
00057 -I if input parameter number I is out of range
00058 1 if answer appears to be lower than lowest
00059 search bound
00060 2 if answer appears to be higher than greatest
00061 search bound
00062 3 if P + Q .ne. 1
00063 4 if X + Y .ne. 1
00064
00065 BOUND <-- Undefined if STATUS is 0
00066
00067 Bound exceeded by parameter number I if STATUS
00068 is negative.
00069
00070 Lower search bound if STATUS is 1.
00071
00072 Upper search bound if STATUS is 2.
00073
00074
00075 Method
00076
00077
00078 Cumulative distribution function (P) is calculated directly by
00079 code associated with the following reference.
00080
00081 DiDinato, A. R. and Morris, A. H. Algorithm 708: Significant
00082 Digit Computation of the Incomplete Beta Function Ratios. ACM
00083 Trans. Math. Softw. 18 (1993), 360-373.
00084
00085 Computation of other parameters involve a seach for a value that
00086 produces the desired value of P. The search relies on the
00087 monotinicity of P with the other parameter.
00088
00089
00090 Note
00091
00092
00093 The beta density is proportional to
00094 t^(A-1) * (1-t)^(B-1)
00095
00096 **********************************************************************/
00097 {
00098 #define tol (1.0e-8)
00099 #define atol (1.0e-50)
00100 #define zero (1.0e-300)
00101 #define inf 1.0e300
00102 #define one 1.0e0
00103 static int K1 = 1;
00104 static double K2 = 0.0e0;
00105 static double K3 = 1.0e0;
00106 static double K8 = 0.5e0;
00107 static double K9 = 5.0e0;
00108 static double fx,xhi,xlo,cum,ccum,xy,pq;
00109 static unsigned long qhi,qleft,qporq;
00110 static double T4,T5,T6,T7,T10,T11,T12,T13,T14,T15;
00111
00112
00113
00114
00115
00116
00117
00118 if(!(*which < 1 || *which > 4)) goto S30;
00119 if(!(*which < 1)) goto S10;
00120 *bound = 1.0e0;
00121 goto S20;
00122 S10:
00123 *bound = 4.0e0;
00124 S20:
00125 *status = -1;
00126 return;
00127 S30:
00128 if(*which == 1) goto S70;
00129
00130
00131
00132 if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
00133 if(!(*p < 0.0e0)) goto S40;
00134 *bound = 0.0e0;
00135 goto S50;
00136 S40:
00137 *bound = 1.0e0;
00138 S50:
00139 *status = -2;
00140 return;
00141 S70:
00142 S60:
00143 if(*which == 1) goto S110;
00144
00145
00146
00147 if(!(*q < 0.0e0 || *q > 1.0e0)) goto S100;
00148 if(!(*q < 0.0e0)) goto S80;
00149 *bound = 0.0e0;
00150 goto S90;
00151 S80:
00152 *bound = 1.0e0;
00153 S90:
00154 *status = -3;
00155 return;
00156 S110:
00157 S100:
00158 if(*which == 2) goto S150;
00159
00160
00161
00162 if(!(*x < 0.0e0 || *x > 1.0e0)) goto S140;
00163 if(!(*x < 0.0e0)) goto S120;
00164 *bound = 0.0e0;
00165 goto S130;
00166 S120:
00167 *bound = 1.0e0;
00168 S130:
00169 *status = -4;
00170 return;
00171 S150:
00172 S140:
00173 if(*which == 2) goto S190;
00174
00175
00176
00177 if(!(*y < 0.0e0 || *y > 1.0e0)) goto S180;
00178 if(!(*y < 0.0e0)) goto S160;
00179 *bound = 0.0e0;
00180 goto S170;
00181 S160:
00182 *bound = 1.0e0;
00183 S170:
00184 *status = -5;
00185 return;
00186 S190:
00187 S180:
00188 if(*which == 3) goto S210;
00189
00190
00191
00192 if(!(*a <= 0.0e0)) goto S200;
00193 *bound = 0.0e0;
00194 *status = -6;
00195 return;
00196 S210:
00197 S200:
00198 if(*which == 4) goto S230;
00199
00200
00201
00202 if(!(*b <= 0.0e0)) goto S220;
00203 *bound = 0.0e0;
00204 *status = -7;
00205 return;
00206 S230:
00207 S220:
00208 if(*which == 1) goto S270;
00209
00210
00211
00212 pq = *p+*q;
00213 if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S260;
00214 if(!(pq < 0.0e0)) goto S240;
00215 *bound = 0.0e0;
00216 goto S250;
00217 S240:
00218 *bound = 1.0e0;
00219 S250:
00220 *status = 3;
00221 return;
00222 S270:
00223 S260:
00224 if(*which == 2) goto S310;
00225
00226
00227
00228 xy = *x+*y;
00229 if(!(fabs(xy-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S300;
00230 if(!(xy < 0.0e0)) goto S280;
00231 *bound = 0.0e0;
00232 goto S290;
00233 S280:
00234 *bound = 1.0e0;
00235 S290:
00236 *status = 4;
00237 return;
00238 S310:
00239 S300:
00240 if(!(*which == 1)) qporq = *p <= *q;
00241
00242
00243
00244
00245 if(1 == *which) {
00246
00247
00248
00249 cumbet(x,y,a,b,p,q);
00250 *status = 0;
00251 }
00252 else if(2 == *which) {
00253
00254
00255
00256 T4 = atol;
00257 T5 = tol;
00258 dstzr(&K2,&K3,&T4,&T5);
00259 if(!qporq) goto S340;
00260 *status = 0;
00261 dzror(status,x,&fx,&xlo,&xhi,&qleft,&qhi);
00262 *y = one-*x;
00263 S320:
00264 if(!(*status == 1)) goto S330;
00265 cumbet(x,y,a,b,&cum,&ccum);
00266 fx = cum-*p;
00267 dzror(status,x,&fx,&xlo,&xhi,&qleft,&qhi);
00268 *y = one-*x;
00269 goto S320;
00270 S330:
00271 goto S370;
00272 S340:
00273 *status = 0;
00274 dzror(status,y,&fx,&xlo,&xhi,&qleft,&qhi);
00275 *x = one-*y;
00276 S350:
00277 if(!(*status == 1)) goto S360;
00278 cumbet(x,y,a,b,&cum,&ccum);
00279 fx = ccum-*q;
00280 dzror(status,y,&fx,&xlo,&xhi,&qleft,&qhi);
00281 *x = one-*y;
00282 goto S350;
00283 S370:
00284 S360:
00285 if(!(*status == -1)) goto S400;
00286 if(!qleft) goto S380;
00287 *status = 1;
00288 *bound = 0.0e0;
00289 goto S390;
00290 S380:
00291 *status = 2;
00292 *bound = 1.0e0;
00293 S400:
00294 S390:
00295 ;
00296 }
00297 else if(3 == *which) {
00298
00299
00300
00301 *a = 5.0e0;
00302 T6 = zero;
00303 T7 = inf;
00304 T10 = atol;
00305 T11 = tol;
00306 dstinv(&T6,&T7,&K8,&K8,&K9,&T10,&T11);
00307 *status = 0;
00308 dinvr(status,a,&fx,&qleft,&qhi);
00309 S410:
00310 if(!(*status == 1)) goto S440;
00311 cumbet(x,y,a,b,&cum,&ccum);
00312 if(!qporq) goto S420;
00313 fx = cum-*p;
00314 goto S430;
00315 S420:
00316 fx = ccum-*q;
00317 S430:
00318 dinvr(status,a,&fx,&qleft,&qhi);
00319 goto S410;
00320 S440:
00321 if(!(*status == -1)) goto S470;
00322 if(!qleft) goto S450;
00323 *status = 1;
00324 *bound = zero;
00325 goto S460;
00326 S450:
00327 *status = 2;
00328 *bound = inf;
00329 S470:
00330 S460:
00331 ;
00332 }
00333 else if(4 == *which) {
00334
00335
00336
00337 *b = 5.0e0;
00338 T12 = zero;
00339 T13 = inf;
00340 T14 = atol;
00341 T15 = tol;
00342 dstinv(&T12,&T13,&K8,&K8,&K9,&T14,&T15);
00343 *status = 0;
00344 dinvr(status,b,&fx,&qleft,&qhi);
00345 S480:
00346 if(!(*status == 1)) goto S510;
00347 cumbet(x,y,a,b,&cum,&ccum);
00348 if(!qporq) goto S490;
00349 fx = cum-*p;
00350 goto S500;
00351 S490:
00352 fx = ccum-*q;
00353 S500:
00354 dinvr(status,b,&fx,&qleft,&qhi);
00355 goto S480;
00356 S510:
00357 if(!(*status == -1)) goto S540;
00358 if(!qleft) goto S520;
00359 *status = 1;
00360 *bound = zero;
00361 goto S530;
00362 S520:
00363 *status = 2;
00364 *bound = inf;
00365 S530:
00366 ;
00367 }
00368 S540:
00369 return;
00370 #undef tol
00371 #undef atol
00372 #undef zero
00373 #undef inf
00374 #undef one
00375 }
|
|
