|
Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
cdf_16.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 | cdfchi (int *which, double *p, double *q, double *x, double *df, int *status, double *bound) |
Define Documentation
Function Documentation
void cdfchi |
( |
int * |
which, |
|
|
double * |
p, |
|
|
double * |
q, |
|
|
double * |
x, |
|
|
double * |
df, |
|
|
int * |
status, |
|
|
double * |
bound |
|
) |
|
|
|
Definition at line 2 of file cdf_16.c.
References cumchi(), dinvr(), dstinv(), p, q, and spmpar().
Referenced by chisq_p2t(), and chisq_t2p().
00025 : 1..3
00026 iwhich = 1 : Calculate P and Q from X and DF
00027 iwhich = 2 : Calculate X from P,Q and DF
00028 iwhich = 3 : Calculate DF from P,Q and X
00029
00030 P <--> The integral from 0 to X of the chi-square
00031 distribution.
00032 Input range: [0, 1].
00033
00034 Q <--> 1-P.
00035 Input range: (0, 1].
00036 P + Q = 1.0.
00037
00038 X <--> Upper limit of integration of the non-central
00039 chi-square distribution.
00040 Input range: [0, +infinity).
00041 Search range: [0,1E300]
00042
00043 DF <--> Degrees of freedom of the
00044 chi-square distribution.
00045 Input range: (0, +infinity).
00046 Search range: [ 1E-300, 1E300]
00047
00048 STATUS <-- 0 if calculation completed correctly
00049 -I if input parameter number I is out of range
00050 1 if answer appears to be lower than lowest
00051 search bound
00052 2 if answer appears to be higher than greatest
00053 search bound
00054 3 if P + Q .ne. 1
00055 10 indicates error returned from cumgam. See
00056 references in cdfgam
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.4.19 of Abramowitz and Stegun, Handbook of
00072 Mathematical Functions (1966) is used to reduce the chisqure
00073 distribution to the incomplete distribution.
00074
00075 Computation of other parameters involve a seach for a value that
00076 produces the desired value of P. The search relies on the
00077 monotinicity of P with the other parameter.
00078
00079 **********************************************************************/
00080 {
00081 #define tol (1.0e-8)
00082 #define atol (1.0e-50)
00083 #define zero (1.0e-300)
00084 #define inf 1.0e300
00085 static int K1 = 1;
00086 static double K2 = 0.0e0;
00087 static double K4 = 0.5e0;
00088 static double K5 = 5.0e0;
00089 static double fx,cum,ccum,pq,porq;
00090 static unsigned long qhi,qleft,qporq;
00091 static double T3,T6,T7,T8,T9,T10,T11;
00092
00093
00094
00095
00096
00097
00098
00099 if(!(*which < 1 || *which > 3)) goto S30;
00100 if(!(*which < 1)) goto S10;
00101 *bound = 1.0e0;
00102 goto S20;
00103 S10:
00104 *bound = 3.0e0;
00105 S20:
00106 *status = -1;
00107 return;
00108 S30:
00109 if(*which == 1) goto S70;
00110
00111
00112
00113 if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
00114 if(!(*p < 0.0e0)) goto S40;
00115 *bound = 0.0e0;
00116 goto S50;
00117 S40:
00118 *bound = 1.0e0;
00119 S50:
00120 *status = -2;
00121 return;
00122 S70:
00123 S60:
00124 if(*which == 1) goto S110;
00125
00126
00127
00128 if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
00129 if(!(*q <= 0.0e0)) goto S80;
00130 *bound = 0.0e0;
00131 goto S90;
00132 S80:
00133 *bound = 1.0e0;
00134 S90:
00135 *status = -3;
00136 return;
00137 S110:
00138 S100:
00139 if(*which == 2) goto S130;
00140
00141
00142
00143 if(!(*x < 0.0e0)) goto S120;
00144 *bound = 0.0e0;
00145 *status = -4;
00146 return;
00147 S130:
00148 S120:
00149 if(*which == 3) goto S150;
00150
00151
00152
00153 if(!(*df <= 0.0e0)) goto S140;
00154 *bound = 0.0e0;
00155 *status = -5;
00156 return;
00157 S150:
00158 S140:
00159 if(*which == 1) goto S190;
00160
00161
00162
00163 pq = *p+*q;
00164 if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S180;
00165 if(!(pq < 0.0e0)) goto S160;
00166 *bound = 0.0e0;
00167 goto S170;
00168 S160:
00169 *bound = 1.0e0;
00170 S170:
00171 *status = 3;
00172 return;
00173 S190:
00174 S180:
00175 if(*which == 1) goto S220;
00176
00177
00178
00179 qporq = *p <= *q;
00180 if(!qporq) goto S200;
00181 porq = *p;
00182 goto S210;
00183 S200:
00184 porq = *q;
00185 S220:
00186 S210:
00187
00188
00189
00190 if(1 == *which) {
00191
00192
00193
00194 *status = 0;
00195 cumchi(x,df,p,q);
00196 if(porq > 1.5e0) {
00197 *status = 10;
00198 return;
00199 }
00200 }
00201 else if(2 == *which) {
00202
00203
00204
00205 *x = 5.0e0;
00206 T3 = inf;
00207 T6 = atol;
00208 T7 = tol;
00209 dstinv(&K2,&T3,&K4,&K4,&K5,&T6,&T7);
00210 *status = 0;
00211 dinvr(status,x,&fx,&qleft,&qhi);
00212 S230:
00213 if(!(*status == 1)) goto S270;
00214 cumchi(x,df,&cum,&ccum);
00215 if(!qporq) goto S240;
00216 fx = cum-*p;
00217 goto S250;
00218 S240:
00219 fx = ccum-*q;
00220 S250:
00221 if(!(fx+porq > 1.5e0)) goto S260;
00222 *status = 10;
00223 return;
00224 S260:
00225 dinvr(status,x,&fx,&qleft,&qhi);
00226 goto S230;
00227 S270:
00228 if(!(*status == -1)) goto S300;
00229 if(!qleft) goto S280;
00230 *status = 1;
00231 *bound = 0.0e0;
00232 goto S290;
00233 S280:
00234 *status = 2;
00235 *bound = inf;
00236 S300:
00237 S290:
00238 ;
00239 }
00240 else if(3 == *which) {
00241
00242
00243
00244 *df = 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,df,&fx,&qleft,&qhi);
00252 S310:
00253 if(!(*status == 1)) goto S350;
00254 cumchi(x,df,&cum,&ccum);
00255 if(!qporq) goto S320;
00256 fx = cum-*p;
00257 goto S330;
00258 S320:
00259 fx = ccum-*q;
00260 S330:
00261 if(!(fx+porq > 1.5e0)) goto S340;
00262 *status = 10;
00263 return;
00264 S340:
00265 dinvr(status,df,&fx,&qleft,&qhi);
00266 goto S310;
00267 S350:
00268 if(!(*status == -1)) goto S380;
00269 if(!qleft) goto S360;
00270 *status = 1;
00271 *bound = zero;
00272 goto S370;
00273 S360:
00274 *status = 2;
00275 *bound = inf;
00276 S370:
00277 ;
00278 }
00279 S380:
00280 return;
00281 #undef tol
00282 #undef atol
00283 #undef zero
00284 #undef inf
00285 }
|
|