00001 #include "cdflib.h"
00002 void cdft(int *which,double *p,double *q,double *t,double *df,
00003           int *status,double *bound)
00004 /**********************************************************************
00005 
00006       void cdft(int *which,double *p,double *q,double *t,double *df,
00007           int *status,double *bound)
00008 
00009                Cumulative Distribution Function
00010                          T distribution
00011 
00012 
00013                               Function
00014 
00015 
00016      Calculates any one parameter of the t distribution given
00017      values for the others.
00018 
00019 
00020                               Arguments
00021 
00022 
00023      WHICH --> Integer indicating which  argument
00024                values is to be calculated from the others.
00025                Legal range: 1..3
00026                iwhich = 1 : Calculate P and Q from T and DF
00027                iwhich = 2 : Calculate T from P,Q and DF
00028                iwhich = 3 : Calculate DF from P,Q and T
00029 
00030         P <--> The integral from -infinity to t of the t-density.
00031                Input range: (0,1].
00032 
00033         Q <--> 1-P.
00034                Input range: (0, 1].
00035                P + Q = 1.0.
00036 
00037         T <--> Upper limit of integration of the t-density.
00038                Input range: ( -infinity, +infinity).
00039                Search range: [ -1E300, 1E300 ]
00040 
00041         DF <--> Degrees of freedom of the t-distribution.
00042                 Input range: (0 , +infinity).
00043                 Search range: [1e-300, 1E10]
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.5.27  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 an incomplete
00069      beta.
00070 
00071      Computation of other parameters involve a seach for a value that
00072      produces  the desired  value  of P.   The search relies  on  the
00073      monotinicity of P with the other parameter.
00074 
00075 **********************************************************************/
00076 {
00077 #define tol (1.0e-8)
00078 #define atol (1.0e-50)
00079 #define zero (1.0e-300)
00080 #define inf 1.0e300
00081 #define maxdf 1.0e10
00082 static int K1 = 1;
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 T2,T3,T6,T7,T8,T9,T10,T11;
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 == 3) goto S130;
00136 /*
00137      DF
00138 */
00139     if(!(*df <= 0.0e0)) goto S120;
00140     *bound = 0.0e0;
00141     *status = -5;
00142     return;
00143 S130:
00144 S120:
00145     if(*which == 1) goto S170;
00146 /*
00147      P + Q
00148 */
00149     pq = *p+*q;
00150     if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S160;
00151     if(!(pq < 0.0e0)) goto S140;
00152     *bound = 0.0e0;
00153     goto S150;
00154 S140:
00155     *bound = 1.0e0;
00156 S150:
00157     *status = 3;
00158     return;
00159 S170:
00160 S160:
00161     if(!(*which == 1)) qporq = *p <= *q;
00162 /*
00163      Select the minimum of P or Q
00164      Calculate ANSWERS
00165 */
00166     if(1 == *which) {
00167 /*
00168      Computing P and Q
00169 */
00170         cumt(t,df,p,q);
00171         *status = 0;
00172     }
00173     else if(2 == *which) {
00174 /*
00175      Computing T
00176      .. Get initial approximation for T
00177 */
00178         *t = dt1(p,q,df);
00179         T2 = -inf;
00180         T3 = inf;
00181         T6 = atol;
00182         T7 = tol;
00183         dstinv(&T2,&T3,&K4,&K4,&K5,&T6,&T7);
00184         *status = 0;
00185         dinvr(status,t,&fx,&qleft,&qhi);
00186 S180:
00187         if(!(*status == 1)) goto S210;
00188         cumt(t,df,&cum,&ccum);
00189         if(!qporq) goto S190;
00190         fx = cum-*p;
00191         goto S200;
00192 S190:
00193         fx = ccum-*q;
00194 S200:
00195         dinvr(status,t,&fx,&qleft,&qhi);
00196         goto S180;
00197 S210:
00198         if(!(*status == -1)) goto S240;
00199         if(!qleft) goto S220;
00200         *status = 1;
00201         *bound = -inf;
00202         goto S230;
00203 S220:
00204         *status = 2;
00205         *bound = inf;
00206 S240:
00207 S230:
00208         ;
00209     }
00210     else if(3 == *which) {
00211 /*
00212      Computing DF
00213 */
00214         *df = 5.0e0;
00215         T8 = zero;
00216         T9 = maxdf;
00217         T10 = atol;
00218         T11 = tol;
00219         dstinv(&T8,&T9,&K4,&K4,&K5,&T10,&T11);
00220         *status = 0;
00221         dinvr(status,df,&fx,&qleft,&qhi);
00222 S250:
00223         if(!(*status == 1)) goto S280;
00224         cumt(t,df,&cum,&ccum);
00225         if(!qporq) goto S260;
00226         fx = cum-*p;
00227         goto S270;
00228 S260:
00229         fx = ccum-*q;
00230 S270:
00231         dinvr(status,df,&fx,&qleft,&qhi);
00232         goto S250;
00233 S280:
00234         if(!(*status == -1)) goto S310;
00235         if(!qleft) goto S290;
00236         *status = 1;
00237         *bound = zero;
00238         goto S300;
00239 S290:
00240         *status = 2;
00241         *bound = maxdf;
00242 S300:
00243         ;
00244     }
00245 S310:
00246     return;
00247 #undef tol
00248 #undef atol
00249 #undef zero
00250 #undef inf
00251 #undef maxdf
00252 } /* END */