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

Diff of /gsl/cdf/tdistinv.c

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

revision 1.1.2.2 by bjg, Mon Jul 14 14:33:22 2003 UTC revision 1.1.2.3 by bjg, Tue Jul 15 17:56:47 2003 UTC
# Line 26  Line 26 
26    
27  #include <stdio.h>  #include <stdio.h>
28    
29    static double
30    inv_cornish_fisher (double z, double nu)
31    {
32      double a = 1 / (nu - 0.5);
33      double b = 48.0 / (a * a);
34    
35      double cf1 = z * (3 + z * z);
36      double cf2 = z * (945 + z * z * (360 + z * z * (63 + z * z * 4)));
37    
38      double y = z - cf1 / b + cf2 / (10 * b * b);
39    
40      double t = GSL_SIGN (z) * sqrt (nu * expm1 (a * y * y));
41    
42      return t;
43    }
44    
45    
46  double  double
47  gsl_cdf_tdist_Pinv (const double p, const double nu)  gsl_cdf_tdist_Pinv (const double p, const double nu)
48  {  {
49    double x;    double x;
50      
51    if (p == 1.0)    if (p == 1.0)
52      {      {
53        return GSL_POSINF;        return GSL_POSINF;
# Line 40  gsl_cdf_tdist_Pinv (const double p, cons Line 57  gsl_cdf_tdist_Pinv (const double p, cons
57        return GSL_NEGINF;        return GSL_NEGINF;
58      }      }
59    
60    if (p < 0.05)  
61      if (p < 0.05)
62      {      {
63                double beta = gsl_sf_beta (0.5, nu / 2);
64          x = -sqrt (nu) * pow (beta * nu * p, -1.0 / nu);
65      }      }
66    else if (p > 0.95)    else if (p > 0.95)
67      {      {
68                double beta = gsl_sf_beta (0.5, nu / 2);
69          x = sqrt (nu) * pow (beta * nu * (1 - p), -1.0 / nu);
70      }      }
71    else    else
72      {      {
73        double xg = gsl_cdf_ugaussian_Pinv (p);        double xg = gsl_cdf_ugaussian_Pinv (p);
74        double x0 = xg + xg*(xg*xg+1.0)/(4.0*nu) ;        x = inv_cornish_fisher (xg, nu);
       x = x0;  
75      }      }
76    
77    
78    {    {
79      double dp, phi;      double dp, phi;
80    
# Line 66  gsl_cdf_tdist_Pinv (const double p, cons Line 86  gsl_cdf_tdist_Pinv (const double p, cons
86        goto end;        goto end;
87    
88      {      {
89        double step = dp / phi;        double lambda = dp / phi;
90          double step0 = lambda;
91          double step1 = ((nu + 1) * x / (x * x + nu)) * (lambda * lambda / 4.0);
92    
93          double step = step0;
94    
95          if (fabs (step1) < fabs (step0))
96            step += step1;
97    
98          printf ("% .18e % .18e % .18e\n", x, step, p - gsl_cdf_tdist_P (x, nu));
99    
       printf("x=%.18e step=%.18e\n", x, step);  
100        x += step;        x += step;
101          
102        if (fabs (step) > 1e-10 * fabs(x))        if (fabs (step) > 1e-10 * fabs (x))
103          goto start;          goto start;
104      }      }
105    }    }
106    
# Line 81  end: Line 109  end:
109    return x;    return x;
110  }  }
111    
112  #if 0  double
113  /*  gsl_cdf_tdist_Qinv (const double p, const double nu)
  * Invert the T distribution. Uses a method as shown in  
  * Statistical Computing, 5.4.2. This method uses an initial  
  * approximation based on Hill's method above, then improves  
  * the initial method with a Taylor series or a Cornish-Fisher  
  * expansion.  
  */  
 static double  
 inv_cornish_fisher (double z, double n)  
 {  
   double ret_val = 0.0;  
   double b;  
   double zsq;  
   double tmp;  
   double c;  
   double u;  
   double d;  
   double zcube;  
   double zfour;  
   double zfive;  
   double zseven;  
   int i;  
   
   tmp = n - 0.5;  
   b = 48.0 * tmp * tmp;  
   if (n > 5.0)  
     {  
       c =  
         96.36 - 16.0 / tmp - 98.0 / (tmp * tmp) +  
         20700.0 / (tmp * tmp * tmp * b);  
     }  
   else  
     {  
       c = 0.3 * (n - 4.5) * (z + 0.6);  
     }  
   zsq = z * z;  
   zcube = zsq * z;  
   zfour = zcube * z;  
   zfive = zfour * z;  
   zseven = zfour * zcube;  
   d = n * M_SQRTPI * gsl_sf_gamma (tmp) / (2.0 * gsl_sf_gamma (tmp + 0.5));  
   u =  
     10.0 * b * (b + c - 2.0 * z - 7.0 * zsq - 5.0 * zcube + 0.05 * d * zfour);  
   tmp = 4.0 * zseven + 63.0 * zfive + 360.0 * zcube + 945.0 * z;  
   ret_val = z - (zcube + 3.0 * z) / b + tmp / u;  
   return ret_val;  
 }  
   
 static double  
 tdist_pdf (const double x, const double nu)  
114  {  {
   double p;  
   
   double lg1 = gsl_sf_lngamma (nu / 2);  
   double lg2 = gsl_sf_lngamma ((nu + 1) / 2);  
   
   p = ((exp (lg2 - lg1) / sqrt (M_PI * nu))  
        * pow ((1 + x * x / nu), -(nu + 1) / 2));  
   return p;  
 }  
   
   double initial_result;  
   double result;  
   double d;  
   double tmp;  
   double tmp2;  
   double z;  
   double zz;  
   double zzz;  
   double y;  
   double tsqn;  
   double a;  
   double method_test;  
   double psi;  
   double psi_prime;  
   double d_psi_dt;  
   double d_psiprime_dt;  
   double tcdf;  
115    double x;    double x;
   double w;  
   double c2;  
   double c3;  
   double c4;  
116    
117    if (prob < 0.0)    if (p == 0.0)
118      {      {
119        return GSL_EDOM;        return GSL_POSINF;
120      }      }
121    if (prob > 1.0)    else if (p == 1.0)
122      {      {
123        return GSL_EDOM;        return GSL_NEGINF;
124      }      }
125    if (nu < 0.0)  
126      {
127        double xg = gsl_cdf_ugaussian_Qinv (p);
128        x = inv_cornish_fisher (xg, nu);
129      }
130    
131      if (x > nu)
132      {      {
133        return GSL_EDOM;        double beta = gsl_sf_beta (0.5, nu / 2);
134          x = sqrt (nu) * pow (beta * nu * p, -1.0 / nu);
135      }      }
136    if (fabs (prob) < GSL_DBL_EPSILON)    else if (x < -nu)
137      {      {
138        return GSL_POSINF;        double beta = gsl_sf_beta (0.5, nu / 2);
139          x = -sqrt (nu) * pow (beta * nu * (1 - p), -1.0 / nu);
140      }      }
141    if (fabs (1.0 - prob) < GSL_DBL_EPSILON)  
142      {
143        double dp, phi;
144    
145      start:
146        dp = -(p - gsl_cdf_tdist_Q (x, nu));
147        phi = gsl_ran_tdist_pdf (x, nu);
148    
149        if (dp == 0.0)
150          goto end;
151    
152      {      {
153        return GSL_NEGINF;        double lambda = dp / phi;
154      }        double step0 = lambda;
155    printf ("prob is %f\t", prob);        double step1 = ((nu + 1) * x / (x * x + nu)) * (lambda * lambda / 4.0);
156    tmp = nu / 2.0;  
157    d = tmp * M_SQRTPI * gsl_sf_gamma (tmp) / gsl_sf_gamma (tmp + 0.5);        double step = step0;
158    method_test = gsl_max (d * prob, 2.0 / nu);  
159    /*        if (fabs (step1) < fabs (step0))
160     * There are two possible initial approximations.          step += step1;
161     * Which is used depends on prob and nu.  
162     */        printf ("% .18e % .18e % .18e\n", x, step, p - gsl_cdf_tdist_P (x, nu));
163    if (method_test > 0.05)  
164      {        x += step;
165        tmp = prob;  
166        a = nu - 0.5;        if (fabs (step) > 1e-10 * fabs (x))
167        /* gsl_cdf_ugaussian_Pinv(tmp) ? */          goto start;
       x = gsl_cdf_ugaussian_Pinv (tmp);  
       y = inv_cornish_fisher (x, nu);  
       tsqn = -1.0 + exp (a * y * y);  
168      }      }
169    else    }
170      {  
171        z = pow (prob * d, 1 / tmp);  end:
172        zz = z * z;  
173        zzz = z * z * z;    return x;
       tsqn = 1 / z + (nu + 1.0) * (-1.0 + z / (2.0 * (nu + 4.0)) +  
                                    nu * zz / (3.0 * (nu + 2.0) * (nu + 6.0))  
                                    + nu * (nu + 3.0) * (2.0 * nu * nu +  
                                                         9.0 * nu -  
                                                         2.0) * zzz / (8.0 *  
                                                                       (nu +  
                                                                        2.0) *  
                                                                       (nu +  
                                                                        2.0) *  
                                                                       (nu +  
                                                                        4.0) *  
                                                                       (nu +  
                                                                        4.0) *  
                                                                       (nu +  
                                                                        8.0)))  
         / (nu + 2);  
   
     }  
   initial_result = sqrt (tsqn * nu);  
   printf ("initial_result = %f \t", initial_result);  
   tcdf = gsl_cdf_tdist_Q (initial_result, nu);  
   tmp = 0.5 * (tcdf - prob);  
   tmp2 = tdist_pdf (initial_result, nu);  
   w = tmp / tmp2;  
   psi = initial_result * (nu + 1.0) / (nu + initial_result * initial_result);  
   psi_prime = (nu + 1.0) * (nu - initial_result * initial_result) /  
     ((nu + initial_result * initial_result) *  
      (nu + initial_result * initial_result));  
   d_psi_dt = 2 * initial_result * nu * (nu + 1.0) /  
     ((nu + initial_result * initial_result) *  
      (nu + initial_result * initial_result));  
   d_psiprime_dt = -2.0 * initial_result * (nu + 1.0) *  
     (3.0 * nu - initial_result * initial_result) /  
     ((nu + initial_result * initial_result) *  
      (nu + initial_result * initial_result) *  
      (nu + initial_result * initial_result));  
   c2 = psi / 2.0;  
   c3 = (2.0 * psi * psi + psi_prime) / 6.0;  
   c4 = 3.0 * psi * (2.0 * psi * psi + psi_prime) +  
     4.0 * psi * d_psi_dt + d_psiprime_dt;  
   c4 /= 24.0;  
   result = initial_result + w + c2 * w * w +  
     c3 * w * w * w + c4 * w * w * w * w;  
   printf ("result = %f\n", result);  
   return result;  
174  }  }
 #endif  

Legend:
Removed from v.1.1.2.2  
changed lines
  Added in v.1.1.2.3

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