/[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.3 by bjg, Tue Jul 15 17:56:47 2003 UTC revision 1.1.2.4 by bjg, Fri Jul 18 13:01:01 2003 UTC
# Line 46  inv_cornish_fisher (double z, double nu) Line 46  inv_cornish_fisher (double z, double nu)
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, ptail;
50    
51    if (p == 1.0)    if (p == 1.0)
52      {      {
# Line 57  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 (nu == 1.0)
   if (p < 0.05)  
61      {      {
62        double beta = gsl_sf_beta (0.5, nu / 2);        x = tan (M_PI * (p - 0.5));
       x = -sqrt (nu) * pow (beta * nu * p, -1.0 / nu);  
63      }      }
64    else if (p > 0.95)    else if (nu == 2.0)
65      {      {
66        double beta = gsl_sf_beta (0.5, nu / 2);        double a = 2 * p - 1;
67        x = sqrt (nu) * pow (beta * nu * (1 - p), -1.0 / nu);        x = a / sqrt (2 * (1 - a * a));
68      }      }
69    else  
70      ptail = (p < 0.5) ? p : 1 - p;
71    
72      if (sqrt (M_PI * nu / 2) * ptail > pow (0.05, nu / 2))
73      {      {
74        double xg = gsl_cdf_ugaussian_Pinv (p);        double xg = gsl_cdf_ugaussian_Pinv (p);
75        x = inv_cornish_fisher (xg, nu);        x = inv_cornish_fisher (xg, nu);
76      }      }
77      else
78        {
79          /* Use an asymptotic expansion of the tail of integral */
80    
81          double beta = gsl_sf_beta (0.5, nu / 2);
82    
83          if (p < 0.5)
84            {
85              x = -sqrt (nu) * pow (beta * nu * p, -1.0 / nu);
86            }
87          else
88            {
89              x = sqrt (nu) * pow (beta * nu * (1 - p), -1.0 / nu);
90            }
91    
92          /* Correct nu -> nu/(1+nu/x^2) in the leading term to account
93             for higher order terms. This avoids overestimating x, which
94             makes the iteration unstable due to the rapidly decreasing
95             tails of the distribution. */
96    
97          x /= sqrt (1 + nu / (x * x));
98        }
99    
100    {    {
101      double dp, phi;      double dp, phi;
# Line 93  gsl_cdf_tdist_Pinv (const double p, cons Line 115  gsl_cdf_tdist_Pinv (const double p, cons
115        double step = step0;        double step = step0;
116    
117        if (fabs (step1) < fabs (step0))        if (fabs (step1) < fabs (step0))
118          step += step1;          {
119              step += step1;
120        printf ("% .18e % .18e % .18e\n", x, step, p - gsl_cdf_tdist_P (x, nu));          }
121    
122        x += step;        if (p > 0.5 && x + step < 0)
123            x /= 2;
124          else if (p < 0.5 && x + step > 0)
125            x /= 2;
126          else
127            x += step;
128    
129        if (fabs (step) > 1e-10 * fabs (x))        if (fabs (step) > 1e-10 * fabs (x))
130          goto start;          goto start;
# Line 112  end: Line 139  end:
139  double  double
140  gsl_cdf_tdist_Qinv (const double p, const double nu)  gsl_cdf_tdist_Qinv (const double p, const double nu)
141  {  {
142    double x;    double x, ptail;
143    
144    if (p == 0.0)    if (p == 0.0)
145      {      {
# Line 123  gsl_cdf_tdist_Qinv (const double p, cons Line 150  gsl_cdf_tdist_Qinv (const double p, cons
150        return GSL_NEGINF;        return GSL_NEGINF;
151      }      }
152    
153    {    if (nu == 1.0)
154      double xg = gsl_cdf_ugaussian_Qinv (p);      {
155      x = inv_cornish_fisher (xg, nu);        x = tan (M_PI * (0.5 - p));
156    }      }
157      else if (nu == 2.0)
158        {
159          double a = 2 * (1 - p) - 1;
160          x = a / sqrt (2 * (1 - a * a));
161        }
162    
163      ptail = (p < 0.5) ? p : 1 - p;
164    
165    if (x > nu)    if (sqrt (M_PI * nu / 2) * ptail > pow (0.05, nu / 2))
166      {      {
167        double beta = gsl_sf_beta (0.5, nu / 2);        double xg = gsl_cdf_ugaussian_Qinv (p);
168        x = sqrt (nu) * pow (beta * nu * p, -1.0 / nu);        x = inv_cornish_fisher (xg, nu);
169      }      }
170    else if (x < -nu)    else
171      {      {
172          /* Use an asymptotic expansion of the tail of integral */
173    
174        double beta = gsl_sf_beta (0.5, nu / 2);        double beta = gsl_sf_beta (0.5, nu / 2);
175        x = -sqrt (nu) * pow (beta * nu * (1 - p), -1.0 / nu);  
176          if (p < 0.5)
177            {
178              x = sqrt (nu) * pow (beta * nu * p, -1.0 / nu);
179            }
180          else
181            {
182              x = -sqrt (nu) * pow (beta * nu * (1 - p), -1.0 / nu);
183            }
184    
185          /* Correct nu -> nu/(1+nu/x^2) in the leading term to account
186             for higher order terms. This avoids overestimating x, which
187             makes the iteration unstable due to the rapidly decreasing
188             tails of the distribution. */
189    
190          x /= sqrt (1 + nu / (x * x));
191      }      }
192    
193    {    {
# Line 157  gsl_cdf_tdist_Qinv (const double p, cons Line 208  gsl_cdf_tdist_Qinv (const double p, cons
208        double step = step0;        double step = step0;
209    
210        if (fabs (step1) < fabs (step0))        if (fabs (step1) < fabs (step0))
211          step += step1;          {
212              step += step1;
213        printf ("% .18e % .18e % .18e\n", x, step, p - gsl_cdf_tdist_P (x, nu));          }
214    
215        x += step;        if (p < 0.5 && x + step < 0)
216            x /= 2;
217          else if (p > 0.5 && x + step > 0)
218            x /= 2;
219          else
220            x += step;
221    
222        if (fabs (step) > 1e-10 * fabs (x))        if (fabs (step) > 1e-10 * fabs (x))
223          goto start;          goto start;

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

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