/[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.3 by jstover, Wed Feb 5 03:45:42 2003 UTC revision 1.1.2.4 by jstover, Thu Feb 6 03:33:31 2003 UTC
# Line 43  Line 43 
43   * needs my own errors to be fixed.)   * needs my own errors to be fixed.)
44   */   */
45  #include <gsl/gsl_cdf.h>  #include <gsl/gsl_cdf.h>
46    #include <gsl/gsl_sf.h>
47  #include <gsl/gsl_math.h>  #include <gsl/gsl_math.h>
48  #include <gsl/gsl_errno.h>  #include <gsl/gsl_errno.h>
49    
50    #define MAXI 40
51    
52  static double cornish_fisher ( double t, double n )  static double cornish_fisher ( double t, double n )
53  {  {
54    double ret_val = 0.0;    double ret_val = 0.0;
# Line 147  static double tdist_pdf (const double x, Line 150  static double tdist_pdf (const double x,
150    
151  /*  /*
152   * This approximation uses different methods depending on the   * This approximation uses different methods depending on the
153   * degrees of freedom (df) and the argument t. An alternate method   * degrees of freedom (nu) and the argument t. An alternate method
154   * could use the incomplete beta function. I didn't choose that   * could use the incomplete beta function. I didn't choose that
155   * method because Kennedy and Gentle state that it is usually inferior to   * method because Kennedy and Gentle state that it is usually inferior to
156   * the method shown here. But "Statistical Computing" was written   * the method shown here. But "Statistical Computing" was written
# Line 156  static double tdist_pdf (const double x, Line 159  static double tdist_pdf (const double x,
159   * this to be true, the method below can be replaced with that   * this to be true, the method below can be replaced with that
160   * of the incomplete beta function.   * of the incomplete beta function.
161   */   */
162  int gsl_cdf_t_e ( double t, double df, gsl_cdf_result * result,  int gsl_cdf_t_e ( double t, double nu, gsl_cdf_result * result,
163                    gsl_cdf_tail_t tail )                    gsl_cdf_tail_t tail )
164  {  {
165    int rc;    int rc;
# Line 167  int gsl_cdf_t_e ( double t, double df, g Line 170  int gsl_cdf_t_e ( double t, double df, g
170    double y;    double y;
171    double diff;    double diff;
172    double ckp2;    double ckp2;
173    double lg1;    gsl_sf_result lg1;
174    double lg2;    gsl_sf_result lg2;
175    double num;    double num;
176    int k;    int k;
177    double u;    double u;
178    double p;    double p;
179    int i;    int i;
180    int idf = (int) df;    int inu = (int) nu;
181    
182    /*    /*
183     * First approximation is also found in Abramowitz     * First approximation is also found in Abramowitz
184     * and Stegun. This is used only for small values of t     * and Stegun. This is used only for small values of t
185     * and small degrees of freedom.     * and small degrees of freedom.
186     */     */
187    y = t / sqrt(df);    y = t / sqrt(nu);
188    if ( fabs(t) < 4.0 )    if ( fabs(t) < 4.0 )
189      {      {
190        diff = fabs(df-1.0);        diff = fabs(nu-1.0);
191        if ( diff < GSL_DBL_EPSILON )        if ( diff < GSL_DBL_EPSILON )
192          {          {
193            q = M_PI_2 * atan(y);            q = M_PI_2 * atan(y);
194          }          }
195        else if ( (df < 21.0) && (df > 1.0))        else if ( (nu < 21.0) && (nu > 1.0))
196          {          {
197            ckp2 = 1.0;            ckp2 = 1.0;
198            b = 1.0 + t*t/df;            b = 1.0 + t*t/nu;
199                        
200            for ( k = idf-2; k > 1; k-=2)            for ( k = inu-2; k > 1; k-=2)
201              {              {
202                ck = 1 + (ckp2 * (k-1)) / ((double)k*b);                ck = 1 + (ckp2 * (k-1)) / ((double)k*b);
203                ckp2 = ck;                ckp2 = ck;
204              }              }
205            if ( (idf%2) == 0 )            if ( (inu%2) == 0 )
206              {              {
207                q = ck * y / sqrt(b);                q = ck * y / sqrt(b);
208              }              }
# Line 216  int gsl_cdf_t_e ( double t, double df, g Line 219  int gsl_cdf_t_e ( double t, double df, g
219                result->val = (1.0+q)/2;                result->val = (1.0+q)/2;
220              }              }
221          }          }
222        else if ( df >= 21.0 )        else if ( nu >= 21.0 )
223          {          {
224            /*            /*
225             * Use the Cornish-Fisher expansion to find             * Use the Cornish-Fisher expansion to find
226             * a point u such that gsl_cdf_gauss(u) = tcdf(t).             * a point u such that gsl_cdf_gauss(u) = tcdf(t).
227             * Approximate the t cdf with gsl_cdf_gauss.             * Approximate the t cdf with gsl_cdf_gauss.
228             */             */
229            u = cornish_fisher ( t, df );            u = cornish_fisher ( t, nu );
230            q = gsl_cdf_gauss ( u, GSL_CDF_UPPER );            q = gsl_cdf_gauss ( u, GSL_CDF_UPPER );
231            if ( tail == GSL_CDF_UPPER )            if ( tail == GSL_CDF_UPPER )
232              {              {
# Line 237  int gsl_cdf_t_e ( double t, double df, g Line 240  int gsl_cdf_t_e ( double t, double df, g
240        else        else
241          {          {
242            GSL_ERROR_VAL ("degrees of freedom are less than 1",            GSL_ERROR_VAL ("degrees of freedom are less than 1",
243                           GSL_EINVAL, df );                           GSL_EINVAL, nu );
244          }          }
245      }      }
246    else    else
# Line 248  int gsl_cdf_t_e ( double t, double df, g Line 251  int gsl_cdf_t_e ( double t, double df, g
251         * to use two different series expansions.         * to use two different series expansions.
252         */         */
253    
254        lg1 = gsl_sf_lngamma (df / 2);  
255        lg2 = gsl_sf_lngamma ((df + 1) / 2);        y = 1/sqrt(1+t*t/nu);
       y = 1/sqrt(1+t*t/df);  
256        num = y;        num = y;
257        p = 0.0;        p = 0.0;
258        for ( i = 2; i < 20; i+=2)  
259          for ( i = 2; (i < MAXI)&&(diff>GSL_DBL_EPSILON); i+=2)
260          {          {
261              diff = p;
262            num *= ((i-1)/i);            num *= ((i-1)/i);
263            p += num / (df+i);            p += num / (nu+i);
264              diff = p - diff;
265            }
266          p = p + 1/nu;
267          rc = gsl_sf_lngamma_e ((nu / 2.0), &lg1);
268          if ( rc == GSL_SUCCESS )
269            {
270              rc = gsl_sf_lngamma_e (((nu + 1) / 2.0), &lg2);
271              if ( rc != GSL_SUCCESS )
272                {
273                  return rc;
274                }
275            }
276          else
277            {
278              return rc;
279          }          }
280        p += 1/df;        diff = (lg2.val) - (lg1.val);
281        p *= pow(y,idf) * exp (lg2 - lg1) / sqrt (M_PI);        p = p * pow(y,nu) * exp (diff) / sqrt (M_PI);
282    
283        if ( tail == GSL_CDF_UPPER )        if ( tail == GSL_CDF_UPPER )
284          {          {
285            result->val = p;            result->val = p;
286          }          }
287        else        else if ( tail == GSL_CDF_LOWER )
288          {          {
289            result->val = 1.0 - p;            result->val = 1.0 - p;
290          }          }

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