/[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.4 by bjg, Fri Jul 18 13:01:01 2003 UTC revision 1.1.2.5 by bjg, Sun Jul 20 15:02:50 2003 UTC
# Line 44  inv_cornish_fisher (double z, double nu) Line 44  inv_cornish_fisher (double z, double nu)
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, ptail;    double x, ptail;
50    
51    if (p == 1.0)    if (P == 1.0)
52      {      {
53        return GSL_POSINF;        return GSL_POSINF;
54      }      }
55    else if (p == 0.0)    else if (P == 0.0)
56      {      {
57        return GSL_NEGINF;        return GSL_NEGINF;
58      }      }
59    
60    if (nu == 1.0)    if (nu == 1.0)
61      {      {
62        x = tan (M_PI * (p - 0.5));        x = tan (M_PI * (P - 0.5));
63      }      }
64    else if (nu == 2.0)    else if (nu == 2.0)
65      {      {
66        double a = 2 * p - 1;        double a = 2 * P - 1;
67        x = a / sqrt (2 * (1 - a * a));        x = a / sqrt (2 * (1 - a * a));
68      }      }
69    
70    ptail = (p < 0.5) ? p : 1 - p;    ptail = (P < 0.5) ? P : 1 - P;
71    
72    if (sqrt (M_PI * nu / 2) * ptail > pow (0.05, nu / 2))    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    else
# Line 80  gsl_cdf_tdist_Pinv (const double p, cons Line 80  gsl_cdf_tdist_Pinv (const double p, cons
80    
81        double beta = gsl_sf_beta (0.5, nu / 2);        double beta = gsl_sf_beta (0.5, nu / 2);
82    
83        if (p < 0.5)        if (P < 0.5)
84          {          {
85            x = -sqrt (nu) * pow (beta * nu * p, -1.0 / nu);            x = -sqrt (nu) * pow (beta * nu * P, -1.0 / nu);
86          }          }
87        else        else
88          {          {
89            x = sqrt (nu) * pow (beta * nu * (1 - p), -1.0 / nu);            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        /* Correct nu -> nu/(1+nu/x^2) in the leading term to account
# Line 98  gsl_cdf_tdist_Pinv (const double p, cons Line 98  gsl_cdf_tdist_Pinv (const double p, cons
98      }      }
99    
100    {    {
101      double dp, phi;      double dP, phi;
102    
103    start:    start:
104      dp = p - gsl_cdf_tdist_P (x, nu);      dP = P - gsl_cdf_tdist_P (x, nu);
105      phi = gsl_ran_tdist_pdf (x, nu);      phi = gsl_ran_tdist_pdf (x, nu);
106    
107      if (dp == 0.0)      if (dP == 0.0)
108        goto end;        goto end;
109    
110      {      {
111        double lambda = dp / phi;        double lambda = dP / phi;
112        double step0 = lambda;        double step0 = lambda;
113        double step1 = ((nu + 1) * x / (x * x + nu)) * (lambda * lambda / 4.0);        double step1 = ((nu + 1) * x / (x * x + nu)) * (lambda * lambda / 4.0);
114    
# Line 119  gsl_cdf_tdist_Pinv (const double p, cons Line 119  gsl_cdf_tdist_Pinv (const double p, cons
119            step += step1;            step += step1;
120          }          }
121    
122        if (p > 0.5 && x + step < 0)        if (P > 0.5 && x + step < 0)
123          x /= 2;          x /= 2;
124        else if (p < 0.5 && x + step > 0)        else if (P < 0.5 && x + step > 0)
125          x /= 2;          x /= 2;
126        else        else
127          x += step;          x += step;
# Line 137  end: Line 137  end:
137  }  }
138    
139  double  double
140  gsl_cdf_tdist_Qinv (const double p, const double nu)  gsl_cdf_tdist_Qinv (const double Q, const double nu)
141  {  {
142    double x, ptail;    double x, qtail;
143    
144    if (p == 0.0)    if (Q == 0.0)
145      {      {
146        return GSL_POSINF;        return GSL_POSINF;
147      }      }
148    else if (p == 1.0)    else if (Q == 1.0)
149      {      {
150        return GSL_NEGINF;        return GSL_NEGINF;
151      }      }
152    
153    if (nu == 1.0)    if (nu == 1.0)
154      {      {
155        x = tan (M_PI * (0.5 - p));        x = tan (M_PI * (0.5 - Q));
156      }      }
157    else if (nu == 2.0)    else if (nu == 2.0)
158      {      {
159        double a = 2 * (1 - p) - 1;        double a = 2 * (1 - Q) - 1;
160        x = a / sqrt (2 * (1 - a * a));        x = a / sqrt (2 * (1 - a * a));
161      }      }
162    
163    ptail = (p < 0.5) ? p : 1 - p;    qtail = (Q < 0.5) ? Q : 1 - Q;
164    
165    if (sqrt (M_PI * nu / 2) * ptail > pow (0.05, nu / 2))    if (sqrt (M_PI * nu / 2) * qtail > pow (0.05, nu / 2))
166      {      {
167        double xg = gsl_cdf_ugaussian_Qinv (p);        double xg = gsl_cdf_ugaussian_Qinv (Q);
168        x = inv_cornish_fisher (xg, nu);        x = inv_cornish_fisher (xg, nu);
169      }      }
170    else    else
# Line 173  gsl_cdf_tdist_Qinv (const double p, cons Line 173  gsl_cdf_tdist_Qinv (const double p, cons
173    
174        double beta = gsl_sf_beta (0.5, nu / 2);        double beta = gsl_sf_beta (0.5, nu / 2);
175    
176        if (p < 0.5)        if (Q < 0.5)
177          {          {
178            x = sqrt (nu) * pow (beta * nu * p, -1.0 / nu);            x = sqrt (nu) * pow (beta * nu * Q, -1.0 / nu);
179          }          }
180        else        else
181          {          {
182            x = -sqrt (nu) * pow (beta * nu * (1 - p), -1.0 / nu);            x = -sqrt (nu) * pow (beta * nu * (1 - Q), -1.0 / nu);
183          }          }
184    
185        /* Correct nu -> nu/(1+nu/x^2) in the leading term to account        /* Correct nu -> nu/(1+nu/x^2) in the leading term to account
# Line 191  gsl_cdf_tdist_Qinv (const double p, cons Line 191  gsl_cdf_tdist_Qinv (const double p, cons
191      }      }
192    
193    {    {
194      double dp, phi;      double dQ, phi;
195    
196    start:    start:
197      dp = -(p - gsl_cdf_tdist_Q (x, nu));      dQ = Q - gsl_cdf_tdist_Q (x, nu);
198      phi = gsl_ran_tdist_pdf (x, nu);      phi = gsl_ran_tdist_pdf (x, nu);
199    
200      if (dp == 0.0)      if (dQ == 0.0)
201        goto end;        goto end;
202    
203      {      {
204        double lambda = dp / phi;        double lambda = - dQ / phi;
205        double step0 = lambda;        double step0 = lambda;
206        double step1 = ((nu + 1) * x / (x * x + nu)) * (lambda * lambda / 4.0);        double step1 = ((nu + 1) * x / (x * x + nu)) * (lambda * lambda / 4.0);
207    
# Line 212  gsl_cdf_tdist_Qinv (const double p, cons Line 212  gsl_cdf_tdist_Qinv (const double p, cons
212            step += step1;            step += step1;
213          }          }
214    
215        if (p < 0.5 && x + step < 0)        if (Q < 0.5 && x + step < 0)
216          x /= 2;          x /= 2;
217        else if (p > 0.5 && x + step > 0)        else if (Q > 0.5 && x + step > 0)
218          x /= 2;          x /= 2;
219        else        else
220          x += step;          x += step;

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

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