/[gsl]/gsl/cdf/t.c
ViewVC logotype

Diff of /gsl/cdf/t.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1.2.13 by jstover, Sun Jun 1 13:58:27 2003 UTC revision 1.1.2.14 by bjg, Fri Jun 20 09:19:45 2003 UTC
# Line 44  Line 44 
44    
45  #define MAXI 40  #define MAXI 40
46    
47  static double cornish_fisher ( double t, double n )  static double
48    cornish_fisher (double t, double n)
49  {  {
50    double ret_val = 0.0;    double ret_val = 0.0;
51    double b;    double b;
52    double z,zsq;    double z, zsq;
53    double tmp;    double tmp;
54    double p[6];    double p[6];
55    int i;    int i;
# Line 63  static double cornish_fisher ( double t, Line 64  static double cornish_fisher ( double t,
64      4960870.65,      4960870.65,
65      37978595.55,      37978595.55,
66      201505390.875,      201505390.875,
67      622437908.625};      622437908.625
68      };
69    const double coeffs5[8] = {    const double coeffs5[8] = {
70      0.2742857142857142857142,      0.2742857142857142857142,
71      4.499047619047619047619,      4.499047619047619047619,
# Line 72  static double cornish_fisher ( double t, Line 74  static double cornish_fisher ( double t,
74      12387.6,      12387.6,
75      101024.55,      101024.55,
76      559494.0,      559494.0,
77      1764959.625};      1764959.625
78    const double coeffs4[6] ={    };
79      const double coeffs4[6] = {
80      0.3047619047619047619048,      0.3047619047619047619048,
81      3.752380952380952380952,      3.752380952380952380952,
82      46.67142857142857142857,      46.67142857142857142857,
83      427.5,      427.5,
84      2587.5,      2587.5,
85      8518.5};      8518.5
86      };
87    const double coeffs3[4] = {    const double coeffs3[4] = {
88      .4,      .4,
89      3.3,      3.3,
90      24.0,      24.0,
91      85.5};      85.5
92      };
93    
94    tmp = n-0.5;    tmp = n - 0.5;
95    z = tmp*log(1+(t*t/n));    z = tmp * log (1 + (t * t / n));
96    z = sqrt(z);    z = sqrt (z);
97    b = 48.0 * tmp * tmp;    b = 48.0 * tmp * tmp;
98    zsq = z*z;    zsq = z * z;
99    
100      p[5] = z * poly_eval(coeffs6, 9, zsq);
101      p[4] = z * poly_eval(coeffs5, 7, zsq);
102      p[3] = z * poly_eval(coeffs4, 5, zsq);
103      p[2] = z * poly_eval(coeffs3, 5, zsq);
104    
105    p[5] = coeffs6[0] * zsq;    p[5] = coeffs6[0] * zsq;
106    for ( i = 1; i < 9; i++)    for (i = 1; i < 9; i++)
107    {      {
108      p[5] = (p[5] + coeffs6[i])*zsq;        p[5] = (p[5] + coeffs6[i]) * zsq;
109    }      }
110    p[5] = z * (p[5] + coeffs6[9]);    p[5] = z * (p[5] + coeffs6[9]);
111      
112    p[4] = coeffs5[0] * zsq;    p[4] = coeffs5[0] * zsq;
113    for(i = 1; i < 7; i++ )    for (i = 1; i < 7; i++)
114      {      {
115        p[4] = (p[4] + coeffs5[i]) * zsq;        p[4] = (p[4] + coeffs5[i]) * zsq;
116      }      }
117    p[4] = z * ( p[4] + coeffs5[7] );    p[4] = z * (p[4] + coeffs5[7]);
118      
119    p[3] = coeffs4[0] * zsq;    p[3] = coeffs4[0] * zsq;
120    for( i = 1; i < 5; i++)    for (i = 1; i < 5; i++)
121      {      {
122        p[3] = (p[3]+coeffs4[i]) * zsq;        p[3] = (p[3] + coeffs4[i]) * zsq;
123      }      }
124    p[3] = z * (p[3] + coeffs4[5]);    p[3] = z * (p[3] + coeffs4[5]);
125    
126    p[2] = coeffs3[0] * zsq;    p[2] = coeffs3[0] * zsq;
127    for(i = 0; i < 3; i++)    for (i = 0; i < 3; i++)
128      {      {
129        p[2] = (p[2] + coeffs3[i]) * zsq;        p[2] = (p[2] + coeffs3[i]) * zsq;
130      }      }
131    p[2] = z * (p[2] + coeffs3[1]);    p[2] = z * (p[2] + coeffs3[1]);  /* BJG: should this be a [3]? */
132    
133    p[1] = z * (zsq +3.0);    p[1] = z * (zsq + 3.0);
134    p[0] = z;    p[0] = z;
135    
136    ret_val = p[5];    ret_val = p[5];
137    for( i = 4; i > -1; i-- )    for (i = 4; i > -1; i--)
138      {      {
139        ret_val = (ret_val/b)+p[i];        ret_val = (ret_val / b) + p[i];
140      }      }
141    return ret_val;    return ret_val;
142  }  }
143    
144  static double tdist_pdf (const double x, const double nu)  static double
145    tdist_pdf (const double x, const double nu)
146  {  {
147    double p;    double p;
148    
# Line 144  static double tdist_pdf (const double x, Line 155  static double tdist_pdf (const double x,
155  }  }
156    
157  /*  /*
  * This approximation uses different methods depending on the  
  * degrees of freedom (nu) and the argument t. An alternate method  
  * could use the incomplete beta function. I didn't choose that  
  * method because Kennedy and Gentle state that it is usually inferior to  
  * the method shown here. But "Statistical Computing" was written  
  * a long time ago, and a recent method for computing the  
  * incomplete beta function may be superior. If experience shows  
  * this to be true, the method below can be replaced with that  
  * of the incomplete beta function.  
  */  
 /* int gsl_cdf_t_e ( double t, double nu) */  
 /* { */  
 /*   int rc; */  
 /*   int stopval; */  
 /*   double b; */  
 /*   double ck; */  
 /*   double q; */  
 /*   double y; */  
 /*   double diff; */  
 /*   double ckp2; */  
 /*   gsl_sf_result lg1; */  
 /*   gsl_sf_result lg2; */  
 /*   double num; */  
 /*   int k; */  
 /*   double u; */  
 /*   double p; */  
 /*   int i; */  
 /*   int inu = (int) nu; */  
   
 /*   /\* */  
 /*    * First approximation is also found in Abramowitz */  
 /*    * and Stegun. This is used only for small values of t */  
 /*    * and small degrees of freedom. */  
 /*    *\/ */  
 /*   y = t / sqrt(nu); */  
 /*   if ( fabs(t) < 4.0 ) */  
 /*     { */  
 /*       diff = fabs(nu-1.0); */  
 /*       if ( diff < GSL_DBL_EPSILON ) */  
 /*      { */  
 /*        q = M_PI_2 * atan(y); */  
 /*      } */  
 /*       else if ( (nu < 21.0) && (nu > 1.0)) */  
 /*      { */  
 /*        ckp2 = 1.0; */  
 /*        b = 1.0 + t*t/nu; */  
             
 /*        for ( k = inu-2; k > 1; k-=2) */  
 /*          { */  
 /*            ck = 1 + (ckp2 * (k-1)) / ((double)k*b); */  
 /*            ckp2 = ck; */  
 /*          } */  
 /*        if ( (inu%2) == 0 ) */  
 /*          { */  
 /*            q = ck * y / sqrt(b); */  
 /*          } */  
 /*        else */  
 /*          { */  
 /*            q = M_PI_2 * (atan(y) + ck * y / b); */  
 /*          } */  
 /*        if ( ((tail == GSL_CDF_UPPER ) && (t > 0.0))||((tail==GSL_CDF_LOWER)&&(t<=0.0))) */  
 /*          { */  
 /*            result->val = (1.0 - q)/2; */  
 /*          } */  
 /*        else */  
 /*          { */  
 /*            result->val = (1.0+q)/2; */  
 /*          } */  
 /*      } */  
 /*       else if ( nu >= 21.0 ) */  
 /*      { */  
           /*  
            * Use the Cornish-Fisher expansion to find  
            * a point u such that gsl_cdf_gauss(u) = tcdf(t).  
            * Approximate the t cdf with gsl_cdf_gauss.  
            */  
 /*        u = cornish_fisher ( t, nu ); */  
 /*        q = gsl_cdf_gauss_Q ( u ); */  
 /*        if ( ((tail == GSL_CDF_UPPER )&&(t>0.0))||((tail==GSL_CDF_LOWER)&&(t<=0.0))) */  
 /*          { */  
 /*            result->val = q; */  
 /*          } */  
 /*        else */  
 /*          { */  
 /*            result->val = 1.0 - q; */  
 /*          } */  
 /*      } */  
 /*       else */  
 /*      { */  
 /*        GSL_ERROR_VAL ("degrees of freedom are less than 1", */  
 /*                       GSL_EINVAL, nu ); */  
 /*      } */  
 /*     } */  
 /*   return rc; */  
 /* } */  
   
 /*  
158   * First approximation is also found in Abramowitz   * First approximation is also found in Abramowitz
159   * and Stegun. This is used only for small values of t   * and Stegun. This is used only for small values of t
160   * and small degrees of freedom.   * and small degrees of freedom.
161   */   */
162  static double gsl_cdf_t_smallx ( const double t, const double nu )  static double
163    gsl_cdf_t_smallx (const double t, const double nu)
164  {  {
165    double y;    double y;
166    double q;    double q;
# Line 257  static double gsl_cdf_t_smallx ( const d Line 172  static double gsl_cdf_t_smallx ( const d
172    int k;    int k;
173    int inu = (int) nu;    int inu = (int) nu;
174    
175    y = t / sqrt(nu);    y = t / sqrt (nu);
176    diff = fabs(nu-1.0);    diff = fabs (nu - 1.0);
177    if ( diff < GSL_DBL_EPSILON )    if (diff < GSL_DBL_EPSILON)
178      {      {
179        q = M_PI_2 * atan(y);        q = M_PI_2 * atan (y);
180      }      }
181    else if ( (nu < 21.0) && (nu > 1.0))    else if ((nu < 21.0) && (nu > 1.0))
182      {      {
183        ckp2 = 1.0;        ckp2 = 1.0;
184        b = 1.0 + t*t/nu;        b = 1.0 + t * t / nu;
185              
186        for ( k = inu-2; k > 1; k-=2)        for (k = inu - 2; k > 1; k -= 2)
187          {          {
188            ck = 1.0 + (ckp2 * (k-1)) / ((double)k*b);            ck = 1.0 + (ckp2 * (k - 1)) / ((double) k * b);
189            ckp2 = ck;            ckp2 = ck;
190          }          }
191        if ( (inu%2) == 0 )        if ((inu % 2) == 0)
192          {          {
193            q = ck * y / sqrt(b);            q = ck * y / sqrt (b);
194          }          }
195        else        else
196          {          {
197            q = M_PI_2 * (atan(y) + ck * y / b);            q = M_PI_2 * (atan (y) + ck * y / b);
198          }          }
 /*       if ( ((tail == GSL_CDF_UPPER ) && (t > 0.0))||((tail==GSL_CDF_LOWER)&&(t<=0.0))) */  
 /*      { */  
 /*        result->val = (1.0 - q)/2; */  
 /*      } */  
 /*       else */  
 /*      { */  
 /*        result->val = (1.0+q)/2; */  
 /*      } */  
199      }      }
200    return q;    return q;
201  }  }
202    
203  /*  /*
204   * Use the Cornish-Fisher expansion to find   * Use the Cornish-Fisher expansion to find
205   * a point u such that gsl_cdf_gauss(u) = tcdf(t).   * a point u such that gsl_cdf_gauss(u) = tcdf(t).
206   * Approximate the t cdf with gsl_cdf_gauss.   * Approximate the t cdf with gsl_cdf_gauss.
207   */   */
208  static double t_cornish_fisher ( const double x, const double nu )  static double
209    t_cornish_fisher (const double x, const double nu)
210  {  {
211    double y;    double y;
212    double q;    double q;
213    double u;    double u;
214    
215    y = x / sqrt(nu);    y = x / sqrt (nu);
216    u = cornish_fisher ( x, nu );    u = cornish_fisher (x, nu);
217    q = gsl_cdf_ugaussian_Q ( u);    q = gsl_cdf_ugaussian_Q (u);
 /*   if ( ((tail == GSL_CDF_UPPER )&&(t>0.0))||((tail==GSL_CDF_LOWER)&&(t<=0.0))) */  
 /*     { */  
 /*       result->val = q; */  
 /*     } */  
 /*   else */  
 /*     { */  
 /*       result->val = 1.0 - q; */  
 /*     } */  
218    return q;    return q;
219  }  }
220    
221  /*  /*
222   * Series approximation for t > 4.0. This needs to be fixed;   * Series approximation for t > 4.0. This needs to be fixed;
223   * it shouldn't subtract the result from 1.0. A better way is   * it shouldn't subtract the result from 1.0. A better way is
# Line 323  static double t_cornish_fisher ( const d Line 225  static double t_cornish_fisher ( const d
225   * means rummaging through Fisher's paper in Metron, v5, 1926,   * means rummaging through Fisher's paper in Metron, v5, 1926,
226   * "Expansion of Student's integral in powers of n^{-1}."   * "Expansion of Student's integral in powers of n^{-1}."
227   */   */
228  static double normal_approx ( const double x, const double nu )  static double
229    normal_approx (const double x, const double nu)
230  {  {
231    double y;    double y;
232    double num;    double num;
# Line 335  static double normal_approx ( const doub Line 238  static double normal_approx ( const doub
238    gsl_sf_result lg1;    gsl_sf_result lg1;
239    gsl_sf_result lg2;    gsl_sf_result lg2;
240    
241    y = 1/sqrt(1+x*x/nu);    y = 1 / sqrt (1 + x * x / nu);
242    num = 1.0;    num = 1.0;
243    q = 0.0;    q = 0.0;
244    diff = 2*GSL_DBL_EPSILON;    diff = 2 * GSL_DBL_EPSILON;
245    for ( i = 2; (i < MAXI)&&(diff>GSL_DBL_EPSILON); i+=2)    for (i = 2; (i < MAXI) && (diff > GSL_DBL_EPSILON); i += 2)
246      {      {
247        diff = q;        diff = q;
248        num *= y*y*(i-1)/i;        num *= y * y * (i - 1) / i;
249        q += num / (nu+i);        q += num / (nu + i);
250        diff = q - diff;        diff = q - diff;
251      }      }
252    q += 1/nu;    q += 1 / nu;
253    rc = gsl_sf_lngamma_e ((nu / 2.0), &lg1);    rc = gsl_sf_lngamma_e ((nu / 2.0), &lg1);
254    if ( rc == GSL_SUCCESS )    if (rc == GSL_SUCCESS)
255      {      {
256        rc = gsl_sf_lngamma_e (((nu + 1) / 2.0), &lg2);        rc = gsl_sf_lngamma_e (((nu + 1) / 2.0), &lg2);
257        if ( rc != GSL_SUCCESS )        if (rc != GSL_SUCCESS)
258          {          {
259            return rc;            return rc;
260          }          }
261      }      }
262    diff = (lg2.val) - (lg1.val);    diff = (lg2.val) - (lg1.val);
263    q *= pow(y,nu) * exp (diff) / sqrt (M_PI);    q *= pow (y, nu) * exp (diff) / sqrt (M_PI);
264    
 /*   if ( ((tail == GSL_CDF_UPPER )&&(t>0.0))||((tail==GSL_CDF_LOWER)&&(t<=0.0))) */  
 /*     { */  
 /*       result->val = q; */  
 /*     } */  
 /*   else if ( tail == GSL_CDF_LOWER ) */  
 /*     { */  
 /*       result->val = 1.0 - q; */  
 /*     } */  
       
265    return q;    return q;
266  }  }
267    
268  double gsl_cdf_t_P ( const double x, const double nu )  double
269    gsl_cdf_t_P (const double x, const double nu)
270  {  {
271    double val = 0.0;    double val = 0.0;
272    double q;    double q;
273    double absx = fabs(x);    double absx = fabs (x);
274      
275    if ( fabs(x) < 4.0 )    if (fabs (x) < 4.0)
276      {      {
277        q = gsl_cdf_t_smallx (x, nu );        q = gsl_cdf_t_smallx (x, nu);
278        if ( x <= 0.0 )        if (x <= 0.0)
279          {          {
280            val = (1.0 - q) / 2.0 ;            val = (1.0 - q) / 2.0;
281          }          }
282        else        else
283          {          {
284            val  = (1.0 + q ) / 2.0;            val = (1.0 + q) / 2.0;
285          }          }
286      }      }
287    else if ( nu >= 21.0 )    else if (nu >= 21.0)
288      {      {
289        q = t_cornish_fisher ( x, nu );        q = t_cornish_fisher (x, nu);
290        if ( x <= 0.0)        if (x <= 0.0)
291          {          {
292            val = q;            val = q;
293          }          }
# Line 403  double gsl_cdf_t_P ( const double x, con Line 298  double gsl_cdf_t_P ( const double x, con
298      }      }
299    else    else
300      {      {
301        val = normal_approx ( x, nu );        val = normal_approx (x, nu);
302        if ( x > 0.0 )        if (x > 0.0)
303          {          {
304            val = 1.0 - val;            val = 1.0 - val;
305          }          }
306      }      }
307    return val;    return val;
308  }  }
309  double gsl_cdf_t_Q ( const double x, const double nu )  
310    double
311    gsl_cdf_t_Q (const double x, const double nu)
312  {  {
313    double val = 0.0;    double val = 0.0;
314    double q;    double q;
315    double absx = fabs(x);    double absx = fabs (x);
316    
317    if ( fabs(x) < 4.0 )    if (fabs (x) < 4.0)
318      {      {
319        q = gsl_cdf_t_smallx ( x, nu );        q = gsl_cdf_t_smallx (x, nu);
320        if ( x > 0.0 )        if (x > 0.0)
321          {          {
322            val = (1.0 - q) / 2.0 ;            val = (1.0 - q) / 2.0;
323          }          }
324        else        else
325          {          {
326            val  = (1.0 + q ) / 2.0;            val = (1.0 + q) / 2.0;
327          }          }
328      }      }
329    else if ( nu >= 21.0 )    else if (nu >= 21.0)
330      {      {
331        val = t_cornish_fisher ( x, nu );        val = t_cornish_fisher (x, nu);
332        if ( x <= 0.0)        if (x <= 0.0)
333          {          {
334            val = 1.0 - val;            val = 1.0 - val;
335          }          }
336      }      }
337    else    else
338      {      {
339        val = normal_approx ( x, nu );        val = normal_approx (x, nu);
340        if ( x <= 0.0 )        if (x <= 0.0)
341          {          {
342            val = 1.0 - val;            val = 1.0 - val;
343          }          }
344      }      }
345    return val;    return val;
346  }  }
347    
348  /*  /*
349   * Invert the T distribution. Uses a method as shown in   * Invert the T distribution. Uses a method as shown in
350   * Statistical Computing, 5.4.2. This method uses an initial   * Statistical Computing, 5.4.2. This method uses an initial
# Line 454  double gsl_cdf_t_Q ( const double x, con Line 352  double gsl_cdf_t_Q ( const double x, con
352   * the initial method with a Taylor series or a Cornish-Fisher   * the initial method with a Taylor series or a Cornish-Fisher
353   * expansion.   * expansion.
354   */   */
355  static double inv_cornish_fisher ( double z, double n )  static double
356    inv_cornish_fisher (double z, double n)
357  {  {
358    double ret_val = 0.0;    double ret_val = 0.0;
359    double b;    double b;
# Line 469  static double inv_cornish_fisher ( doubl Line 368  static double inv_cornish_fisher ( doubl
368    double zseven;    double zseven;
369    int i;    int i;
370    
371    tmp = n-0.5;    tmp = n - 0.5;
372    b = 48.0 * tmp * tmp;    b = 48.0 * tmp * tmp;
373    if ( n > 5.0 )    if (n > 5.0)
374      {      {
375        c = 96.36 - 16.0 / tmp - 98.0 / (tmp*tmp) + 20700.0 / (tmp*tmp*tmp*b);        c =
376            96.36 - 16.0 / tmp - 98.0 / (tmp * tmp) +
377            20700.0 / (tmp * tmp * tmp * b);
378      }      }
379    else    else
380      {      {
381        c = 0.3 * (n - 4.5) * (z + 0.6);        c = 0.3 * (n - 4.5) * (z + 0.6);
382      }      }
383    zsq = z*z;    zsq = z * z;
384    zcube = zsq*z;    zcube = zsq * z;
385    zfour = zcube*z;    zfour = zcube * z;
386    zfive = zfour*z;    zfive = zfour * z;
387    zseven = zfour*zcube;    zseven = zfour * zcube;
388    d = n * M_SQRTPI * gsl_sf_gamma ( tmp )      d = n * M_SQRTPI * gsl_sf_gamma (tmp) / (2.0 * gsl_sf_gamma (tmp + 0.5));
389      / (2.0*gsl_sf_gamma ( tmp + 0.5 ));    u =
390    u = 10.0 * b * ( b + c - 2.0*z - 7.0 * zsq - 5.0*zcube + 0.05*d*zfour);      10.0 * b * (b + c - 2.0 * z - 7.0 * zsq - 5.0 * zcube + 0.05 * d * zfour);
391    tmp = 4.0*zseven + 63.0*zfive + 360.0*zcube + 945.0*z;    tmp = 4.0 * zseven + 63.0 * zfive + 360.0 * zcube + 945.0 * z;
392    ret_val = z - (zcube+3.0*z)/b + tmp/u;    ret_val = z - (zcube + 3.0 * z) / b + tmp / u;
393    return ret_val;    return ret_val;
394  }  }
395    
396  double gsl_cdf_ut_P_inv (double prob, double nu)  double
397    gsl_cdf_ut_P_inv (double prob, double nu)
398  {  {
399    double initial_result;    double initial_result;
400    double result;    double result;
# Line 517  double gsl_cdf_ut_P_inv (double prob, do Line 419  double gsl_cdf_ut_P_inv (double prob, do
419    double c3;    double c3;
420    double c4;    double c4;
421    
422    if ( prob < 0.0 )    if (prob < 0.0)
423      {      {
424        return GSL_EDOM;        return GSL_EDOM;
425      }      }
426    if ( prob > 1.0 )    if (prob > 1.0)
427      {      {
428        return GSL_EDOM;        return GSL_EDOM;
429      }      }
430    if ( nu < 0.0 )    if (nu < 0.0)
431      {      {
432        return GSL_EDOM;        return GSL_EDOM;
433      }      }
434    if ( fabs(prob) < GSL_DBL_EPSILON )    if (fabs (prob) < GSL_DBL_EPSILON)
435      {      {
436        return GSL_POSINF;        return GSL_POSINF;
437      }      }
438    if ( fabs ( 1.0 - prob ) < GSL_DBL_EPSILON )    if (fabs (1.0 - prob) < GSL_DBL_EPSILON)
439      {      {
440        return GSL_NEGINF;        return GSL_NEGINF;
441      }      }
442    printf("prob is %f\t",prob);    printf ("prob is %f\t", prob);
443    tmp = nu / 2.0;    tmp = nu / 2.0;
444    d = tmp * M_SQRTPI * gsl_sf_gamma ( tmp )      d = tmp * M_SQRTPI * gsl_sf_gamma (tmp) / gsl_sf_gamma (tmp + 0.5);
445      / gsl_sf_gamma ( tmp + 0.5 );    method_test = gsl_max (d * prob, 2.0 / nu);
   method_test = gsl_max ( d*prob, 2.0 / nu);  
446    /*    /*
447     * There are two possible initial approximations.     * There are two possible initial approximations.
448     * Which is used depends on prob and nu.     * Which is used depends on prob and nu.
449     */     */
450    if ( method_test > 0.05 )    if (method_test > 0.05)
451      {      {
452        tmp = prob;        tmp = prob;
453        a = nu - 0.5;        a = nu - 0.5;
454        /* gsl_cdf_ugaussian_P_inv(tmp) ? */        /* gsl_cdf_ugaussian_P_inv(tmp) ? */
455        x = gsl_cdf_ugaussian_P_inv (tmp );        x = gsl_cdf_ugaussian_P_inv (tmp);
456        y = inv_cornish_fisher ( x, nu );        y = inv_cornish_fisher (x, nu);
457        tsqn = -1.0 + exp ( a*y*y);        tsqn = -1.0 + exp (a * y * y);
458      }      }
459    else    else
460      {      {
461        z = pow ( prob * d, 1 / tmp );        z = pow (prob * d, 1 / tmp);
462        zz = z * z;        zz = z * z;
463        zzz = z*z*z;        zzz = z * z * z;
464        tsqn = 1/z + (nu + 1.0)*(-1.0 + z/(2.0 * (nu + 4.0)) +        tsqn = 1 / z + (nu + 1.0) * (-1.0 + z / (2.0 * (nu + 4.0)) +
465                               nu * zz / (3.0 * (nu+2.0)*(nu+6.0))                                     nu * zz / (3.0 * (nu + 2.0) * (nu + 6.0))
466                               + nu*(nu+3.0)*(2.0*nu*nu+9.0*nu-2.0)*zzz                                     + nu * (nu + 3.0) * (2.0 * nu * nu +
467                               /(8.0 * (nu + 2.0) * (nu + 2.0)                                                          9.0 * nu -
468                                 * (nu + 4.0) * (nu + 4.0) *                                                          2.0) * zzz / (8.0 *
469                                 (nu + 8.0)))                                                                        (nu +
470          /(nu + 2);                                                                         2.0) *
471                                                                          (nu +
472      }                                                                         2.0) *
473    initial_result = sqrt(tsqn * nu);                                                                        (nu +
474    printf("initial_result = %f \t",initial_result);                                                                         4.0) *
475    tcdf = gsl_cdf_t_Q( initial_result, nu );                                                                        (nu +
476    tmp = 0.5 * ( tcdf - prob );                                                                         4.0) *
477    tmp2 = tdist_pdf ( initial_result, nu );                                                                        (nu +
478                                                                           8.0)))
479            / (nu + 2);
480    
481        }
482      initial_result = sqrt (tsqn * nu);
483      printf ("initial_result = %f \t", initial_result);
484      tcdf = gsl_cdf_t_Q (initial_result, nu);
485      tmp = 0.5 * (tcdf - prob);
486      tmp2 = tdist_pdf (initial_result, nu);
487    w = tmp / tmp2;    w = tmp / tmp2;
488    psi = initial_result * (nu+1.0) /    psi = initial_result * (nu + 1.0) / (nu + initial_result * initial_result);
489      (nu + initial_result * initial_result);    psi_prime = (nu + 1.0) * (nu - initial_result * initial_result) /
490    psi_prime = (nu+1.0) * (nu - initial_result * initial_result) /      ((nu + initial_result * initial_result) *
     ((nu + initial_result*initial_result) *  
491       (nu + initial_result * initial_result));       (nu + initial_result * initial_result));
492    d_psi_dt = 2 * initial_result * nu * (nu + 1.0) /    d_psi_dt = 2 * initial_result * nu * (nu + 1.0) /
493      ((nu + initial_result * initial_result) *      ((nu + initial_result * initial_result) *
# Line 588  double gsl_cdf_ut_P_inv (double prob, do Line 497  double gsl_cdf_ut_P_inv (double prob, do
497      ((nu + initial_result * initial_result) *      ((nu + initial_result * initial_result) *
498       (nu + initial_result * initial_result) *       (nu + initial_result * initial_result) *
499       (nu + initial_result * initial_result));       (nu + initial_result * initial_result));
500    c2 = psi/2.0;    c2 = psi / 2.0;
501    c3 = (2.0 * psi * psi + psi_prime) / 6.0;    c3 = (2.0 * psi * psi + psi_prime) / 6.0;
502    c4 = 3.0 * psi * (2.0 * psi * psi + psi_prime) +    c4 = 3.0 * psi * (2.0 * psi * psi + psi_prime) +
503      4.0 * psi * d_psi_dt + d_psiprime_dt;      4.0 * psi * d_psi_dt + d_psiprime_dt;
504    c4 /= 24.0;    c4 /= 24.0;
505    result = initial_result + w + c2 * w * w +    result = initial_result + w + c2 * w * w +
506      c3 * w * w * w + c4 * w * w * w * w;      c3 * w * w * w + c4 * w * w * w * w;
507    printf("result = %f\n",result);    printf ("result = %f\n", result);
508    return result;    return result;
509  }  }

Legend:
Removed from v.1.1.2.13  
changed lines
  Added in v.1.1.2.14

savannah-hackers-public@gnu.org
ViewVC Help
Powered by ViewVC 1.1.26