Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/ALGLIB/studenttdistr.cs @ 2601

Last change on this file since 2601 was 2563, checked in by gkronber, 15 years ago

Updated ALGLIB to latest version. #751 (Plugin for for data-modeling with ANN (integrated into CEDMA))

File size: 7.3 KB
Line 
1/*************************************************************************
2Cephes Math Library Release 2.8:  June, 2000
3Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
4
5Contributors:
6    * Sergey Bochkanov (ALGLIB project). Translation from C to
7      pseudocode.
8
9See subroutines comments for additional copyrights.
10
11>>> SOURCE LICENSE >>>
12This program is free software; you can redistribute it and/or modify
13it under the terms of the GNU General Public License as published by
14the Free Software Foundation (www.fsf.org); either version 2 of the
15License, or (at your option) any later version.
16
17This program is distributed in the hope that it will be useful,
18but WITHOUT ANY WARRANTY; without even the implied warranty of
19MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
20GNU General Public License for more details.
21
22A copy of the GNU General Public License is available at
23http://www.fsf.org/licensing/licenses
24
25>>> END OF LICENSE >>>
26*************************************************************************/
27
28using System;
29
30namespace alglib
31{
32    public class studenttdistr
33    {
34        /*************************************************************************
35        Student's t distribution
36
37        Computes the integral from minus infinity to t of the Student
38        t distribution with integer k > 0 degrees of freedom:
39
40                                             t
41                                             -
42                                            | |
43                     -                      |         2   -(k+1)/2
44                    | ( (k+1)/2 )           |  (     x   )
45              ----------------------        |  ( 1 + --- )        dx
46                            -               |  (      k  )
47              sqrt( k pi ) | ( k/2 )        |
48                                          | |
49                                           -
50                                          -inf.
51
52        Relation to incomplete beta integral:
53
54               1 - stdtr(k,t) = 0.5 * incbet( k/2, 1/2, z )
55        where
56               z = k/(k + t**2).
57
58        For t < -2, this is the method of computation.  For higher t,
59        a direct method is derived from integration by parts.
60        Since the function is symmetric about t=0, the area under the
61        right tail of the density is found by calling the function
62        with -t instead of t.
63
64        ACCURACY:
65
66        Tested at random 1 <= k <= 25.  The "domain" refers to t.
67                             Relative error:
68        arithmetic   domain     # trials      peak         rms
69           IEEE     -100,-2      50000       5.9e-15     1.4e-15
70           IEEE     -2,100      500000       2.7e-15     4.9e-17
71
72        Cephes Math Library Release 2.8:  June, 2000
73        Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
74        *************************************************************************/
75        public static double studenttdistribution(int k,
76            double t)
77        {
78            double result = 0;
79            double x = 0;
80            double rk = 0;
81            double z = 0;
82            double f = 0;
83            double tz = 0;
84            double p = 0;
85            double xsqk = 0;
86            int j = 0;
87
88            System.Diagnostics.Debug.Assert(k>0, "Domain error in StudentTDistribution");
89            if( (double)(t)==(double)(0) )
90            {
91                result = 0.5;
92                return result;
93            }
94            if( (double)(t)<(double)(-2.0) )
95            {
96                rk = k;
97                z = rk/(rk+t*t);
98                result = 0.5*ibetaf.incompletebeta(0.5*rk, 0.5, z);
99                return result;
100            }
101            if( (double)(t)<(double)(0) )
102            {
103                x = -t;
104            }
105            else
106            {
107                x = t;
108            }
109            rk = k;
110            z = 1.0+x*x/rk;
111            if( k%2!=0 )
112            {
113                xsqk = x/Math.Sqrt(rk);
114                p = Math.Atan(xsqk);
115                if( k>1 )
116                {
117                    f = 1.0;
118                    tz = 1.0;
119                    j = 3;
120                    while( j<=k-2 & (double)(tz/f)>(double)(AP.Math.MachineEpsilon) )
121                    {
122                        tz = tz*((j-1)/(z*j));
123                        f = f+tz;
124                        j = j+2;
125                    }
126                    p = p+f*xsqk/z;
127                }
128                p = p*2.0/Math.PI;
129            }
130            else
131            {
132                f = 1.0;
133                tz = 1.0;
134                j = 2;
135                while( j<=k-2 & (double)(tz/f)>(double)(AP.Math.MachineEpsilon) )
136                {
137                    tz = tz*((j-1)/(z*j));
138                    f = f+tz;
139                    j = j+2;
140                }
141                p = f*x/Math.Sqrt(z*rk);
142            }
143            if( (double)(t)<(double)(0) )
144            {
145                p = -p;
146            }
147            result = 0.5+0.5*p;
148            return result;
149        }
150
151
152        /*************************************************************************
153        Functional inverse of Student's t distribution
154
155        Given probability p, finds the argument t such that stdtr(k,t)
156        is equal to p.
157
158        ACCURACY:
159
160        Tested at random 1 <= k <= 100.  The "domain" refers to p:
161                             Relative error:
162        arithmetic   domain     # trials      peak         rms
163           IEEE    .001,.999     25000       5.7e-15     8.0e-16
164           IEEE    10^-6,.001    25000       2.0e-12     2.9e-14
165
166        Cephes Math Library Release 2.8:  June, 2000
167        Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
168        *************************************************************************/
169        public static double invstudenttdistribution(int k,
170            double p)
171        {
172            double result = 0;
173            double t = 0;
174            double rk = 0;
175            double z = 0;
176            int rflg = 0;
177
178            System.Diagnostics.Debug.Assert(k>0 & (double)(p)>(double)(0) & (double)(p)<(double)(1), "Domain error in InvStudentTDistribution");
179            rk = k;
180            if( (double)(p)>(double)(0.25) & (double)(p)<(double)(0.75) )
181            {
182                if( (double)(p)==(double)(0.5) )
183                {
184                    result = 0;
185                    return result;
186                }
187                z = 1.0-2.0*p;
188                z = ibetaf.invincompletebeta(0.5, 0.5*rk, Math.Abs(z));
189                t = Math.Sqrt(rk*z/(1.0-z));
190                if( (double)(p)<(double)(0.5) )
191                {
192                    t = -t;
193                }
194                result = t;
195                return result;
196            }
197            rflg = -1;
198            if( (double)(p)>=(double)(0.5) )
199            {
200                p = 1.0-p;
201                rflg = 1;
202            }
203            z = ibetaf.invincompletebeta(0.5*rk, 0.5, 2.0*p);
204            if( (double)(AP.Math.MaxRealNumber*z)<(double)(rk) )
205            {
206                result = rflg*AP.Math.MaxRealNumber;
207                return result;
208            }
209            t = Math.Sqrt(rk/z-rk);
210            result = rflg*t;
211            return result;
212        }
213    }
214}
Note: See TracBrowser for help on using the repository browser.