source: sasview/src/sas/models/c_extension/cephes/stdtr.c @ c81f9b6

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since c81f9b6 was 79492222, checked in by krzywon, 10 years ago

Changed the file and folder names to remove all SANS references.

  • Property mode set to 100644
File size: 4.0 KB
Line 
1/*                                                      stdtr.c
2 *
3 *      Student's t distribution
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * double t, stdtr();
10 * short k;
11 *
12 * y = stdtr( k, t );
13 *
14 *
15 * DESCRIPTION:
16 *
17 * Computes the integral from minus infinity to t of the Student
18 * t distribution with integer k > 0 degrees of freedom:
19 *
20 *                                      t
21 *                                      -
22 *                                     | |
23 *              -                      |         2   -(k+1)/2
24 *             | ( (k+1)/2 )           |  (     x   )
25 *       ----------------------        |  ( 1 + --- )        dx
26 *                     -               |  (      k  )
27 *       sqrt( k pi ) | ( k/2 )        |
28 *                                   | |
29 *                                    -
30 *                                   -inf.
31 *
32 * Relation to incomplete beta integral:
33 *
34 *        1 - stdtr(k,t) = 0.5 * incbet( k/2, 1/2, z )
35 * where
36 *        z = k/(k + t**2).
37 *
38 * For t < -2, this is the method of computation.  For higher t,
39 * a direct method is derived from integration by parts.
40 * Since the function is symmetric about t=0, the area under the
41 * right tail of the density is found by calling the function
42 * with -t instead of t.
43 *
44 * ACCURACY:
45 *
46 * Tested at random 1 <= k <= 25.  The "domain" refers to t.
47 *                      Relative error:
48 * arithmetic   domain     # trials      peak         rms
49 *    IEEE     -100,-2      50000       5.9e-15     1.4e-15
50 *    IEEE     -2,100      500000       2.7e-15     4.9e-17
51 */
52
53/*                                                      stdtri.c
54 *
55 *      Functional inverse of Student's t distribution
56 *
57 *
58 *
59 * SYNOPSIS:
60 *
61 * double p, t, stdtri();
62 * int k;
63 *
64 * t = stdtri( k, p );
65 *
66 *
67 * DESCRIPTION:
68 *
69 * Given probability p, finds the argument t such that stdtr(k,t)
70 * is equal to p.
71 *
72 * ACCURACY:
73 *
74 * Tested at random 1 <= k <= 100.  The "domain" refers to p:
75 *                      Relative error:
76 * arithmetic   domain     # trials      peak         rms
77 *    IEEE    .001,.999     25000       5.7e-15     8.0e-16
78 *    IEEE    10^-6,.001    25000       2.0e-12     2.9e-14
79 */
80
81
82/*
83Cephes Math Library Release 2.8:  June, 2000
84Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
85*/
86
87#include "mconf.h"
88
89extern double PI, MACHEP, MAXNUM;
90#ifdef ANSIPROT
91extern double sqrt ( double );
92extern double atan ( double );
93extern double incbet ( double, double, double );
94extern double incbi ( double, double, double );
95extern double fabs ( double );
96#else
97double sqrt(), atan(), incbet(), incbi(), fabs();
98#endif
99
100double stdtr( k, t )
101int k;
102double t;
103{
104double x, rk, z, f, tz, p, xsqk;
105int j;
106
107if( k <= 0 )
108        {
109        mtherr( "stdtr", DOMAIN );
110        return(0.0);
111        }
112
113if( t == 0 )
114        return( 0.5 );
115
116if( t < -2.0 )
117        {
118        rk = k;
119        z = rk / (rk + t * t);
120        p = 0.5 * incbet( 0.5*rk, 0.5, z );
121        return( p );
122        }
123
124/*      compute integral from -t to + t */
125
126if( t < 0 )
127        x = -t;
128else
129        x = t;
130
131rk = k; /* degrees of freedom */
132z = 1.0 + ( x * x )/rk;
133
134/* test if k is odd or even */
135if( (k & 1) != 0)
136        {
137
138        /*      computation for odd k   */
139
140        xsqk = x/sqrt(rk);
141        p = atan( xsqk );
142        if( k > 1 )
143                {
144                f = 1.0;
145                tz = 1.0;
146                j = 3;
147                while(  (j<=(k-2)) && ( (tz/f) > MACHEP )  )
148                        {
149                        tz *= (j-1)/( z * j );
150                        f += tz;
151                        j += 2;
152                        }
153                p += f * xsqk/z;
154                }
155        p *= 2.0/PI;
156        }
157
158
159else
160        {
161
162        /*      computation for even k  */
163
164        f = 1.0;
165        tz = 1.0;
166        j = 2;
167
168        while(  ( j <= (k-2) ) && ( (tz/f) > MACHEP )  )
169                {
170                tz *= (j - 1)/( z * j );
171                f += tz;
172                j += 2;
173                }
174        p = f * x/sqrt(z*rk);
175        }
176
177/*      common exit     */
178
179
180if( t < 0 )
181        p = -p; /* note destruction of relative accuracy */
182
183        p = 0.5 + 0.5 * p;
184return(p);
185}
186
187double stdtri( k, p )
188int k;
189double p;
190{
191double t, rk, z;
192int rflg;
193
194if( k <= 0 || p <= 0.0 || p >= 1.0 )
195        {
196        mtherr( "stdtri", DOMAIN );
197        return(0.0);
198        }
199
200rk = k;
201
202if( p > 0.25 && p < 0.75 )
203        {
204        if( p == 0.5 )
205                return( 0.0 );
206        z = 1.0 - 2.0 * p;
207        z = incbi( 0.5, 0.5*rk, fabs(z) );
208        t = sqrt( rk*z/(1.0-z) );
209        if( p < 0.5 )
210                t = -t;
211        return( t );
212        }
213rflg = -1;
214if( p >= 0.5)
215        {
216        p = 1.0 - p;
217        rflg = 1;
218        }
219z = incbi( 0.5*rk, 0.5, 2.0*p );
220
221if( MAXNUM * z < rk )
222        return(rflg* MAXNUM);
223t = sqrt( rk/z - rk );
224return( rflg * t );
225}
Note: See TracBrowser for help on using the repository browser.