Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/ALGLIB/normaldistr.cs @ 2474

Last change on this file since 2474 was 2445, checked in by gkronber, 15 years ago

Fixed #787 (LinearRegressionOperator uses leastsquares function of ALGLIB instead of linearregression function)
File size: 13.0 KB
Line 
1/*************************************************************************
2Cephes Math Library Release 2.8:  June, 2000
3Copyright 1984, 1987, 1988, 1992, 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 normaldistr
33    {
34        /*************************************************************************
35        Error function
36
37        The integral is
38
39                                  x
40                                   -
41                        2         | |          2
42          erf(x)  =  --------     |    exp( - t  ) dt.
43                     sqrt(pi)   | |
44                                 -
45                                  0
46
47        For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise
48        erf(x) = 1 - erfc(x).
49
50
51        ACCURACY:
52
53                             Relative error:
54        arithmetic   domain     # trials      peak         rms
55           IEEE      0,1         30000       3.7e-16     1.0e-16
56
57        Cephes Math Library Release 2.8:  June, 2000
58        Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
59        *************************************************************************/
60        public static double erf(double x)
61        {
62            double result = 0;
63            double xsq = 0;
64            double s = 0;
65            double p = 0;
66            double q = 0;
67
68            s = Math.Sign(x);
69            x = Math.Abs(x);
70            if( x<0.5 )
71            {
72                xsq = x*x;
73                p = 0.007547728033418631287834;
74                p = 0.288805137207594084924010+xsq*p;
75                p = 14.3383842191748205576712+xsq*p;
76                p = 38.0140318123903008244444+xsq*p;
77                p = 3017.82788536507577809226+xsq*p;
78                p = 7404.07142710151470082064+xsq*p;
79                p = 80437.3630960840172832162+xsq*p;
80                q = 0.0;
81                q = 1.00000000000000000000000+xsq*q;
82                q = 38.0190713951939403753468+xsq*q;
83                q = 658.070155459240506326937+xsq*q;
84                q = 6379.60017324428279487120+xsq*q;
85                q = 34216.5257924628539769006+xsq*q;
86                q = 80437.3630960840172826266+xsq*q;
87                result = s*1.1283791670955125738961589031*x*p/q;
88                return result;
89            }
90            if( x>=10 )
91            {
92                result = s;
93                return result;
94            }
95            result = s*(1-erfc(x));
96            return result;
97        }
98
99
100        /*************************************************************************
101        Complementary error function
102
103         1 - erf(x) =
104
105                                  inf.
106                                    -
107                         2         | |          2
108          erfc(x)  =  --------     |    exp( - t  ) dt
109                      sqrt(pi)   | |
110                                  -
111                                   x
112
113
114        For small x, erfc(x) = 1 - erf(x); otherwise rational
115        approximations are computed.
116
117
118        ACCURACY:
119
120                             Relative error:
121        arithmetic   domain     # trials      peak         rms
122           IEEE      0,26.6417   30000       5.7e-14     1.5e-14
123
124        Cephes Math Library Release 2.8:  June, 2000
125        Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
126        *************************************************************************/
127        public static double erfc(double x)
128        {
129            double result = 0;
130            double p = 0;
131            double q = 0;
132
133            if( x<0 )
134            {
135                result = 2-erfc(-x);
136                return result;
137            }
138            if( x<0.5 )
139            {
140                result = 1.0-erf(x);
141                return result;
142            }
143            if( x>=10 )
144            {
145                result = 0;
146                return result;
147            }
148            p = 0.0;
149            p = 0.5641877825507397413087057563+x*p;
150            p = 9.675807882987265400604202961+x*p;
151            p = 77.08161730368428609781633646+x*p;
152            p = 368.5196154710010637133875746+x*p;
153            p = 1143.262070703886173606073338+x*p;
154            p = 2320.439590251635247384768711+x*p;
155            p = 2898.0293292167655611275846+x*p;
156            p = 1826.3348842295112592168999+x*p;
157            q = 1.0;
158            q = 17.14980943627607849376131193+x*q;
159            q = 137.1255960500622202878443578+x*q;
160            q = 661.7361207107653469211984771+x*q;
161            q = 2094.384367789539593790281779+x*q;
162            q = 4429.612803883682726711528526+x*q;
163            q = 6089.5424232724435504633068+x*q;
164            q = 4958.82756472114071495438422+x*q;
165            q = 1826.3348842295112595576438+x*q;
166            result = Math.Exp(-AP.Math.Sqr(x))*p/q;
167            return result;
168        }
169
170
171        /*************************************************************************
172        Normal distribution function
173
174        Returns the area under the Gaussian probability density
175        function, integrated from minus infinity to x:
176
177                                   x
178                                    -
179                          1        | |          2
180           ndtr(x)  = ---------    |    exp( - t /2 ) dt
181                      sqrt(2pi)  | |
182                                  -
183                                 -inf.
184
185                    =  ( 1 + erf(z) ) / 2
186                    =  erfc(z) / 2
187
188        where z = x/sqrt(2). Computation is via the functions
189        erf and erfc.
190
191
192        ACCURACY:
193
194                             Relative error:
195        arithmetic   domain     # trials      peak         rms
196           IEEE     -13,0        30000       3.4e-14     6.7e-15
197
198        Cephes Math Library Release 2.8:  June, 2000
199        Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
200        *************************************************************************/
201        public static double normaldistribution(double x)
202        {
203            double result = 0;
204
205            result = 0.5*(erf(x/1.41421356237309504880)+1);
206            return result;
207        }
208
209
210        /*************************************************************************
211        Inverse of the error function
212
213        Cephes Math Library Release 2.8:  June, 2000
214        Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
215        *************************************************************************/
216        public static double inverf(double e)
217        {
218            double result = 0;
219
220            result = invnormaldistribution(0.5*(e+1))/Math.Sqrt(2);
221            return result;
222        }
223
224
225        /*************************************************************************
226        Inverse of Normal distribution function
227
228        Returns the argument, x, for which the area under the
229        Gaussian probability density function (integrated from
230        minus infinity to x) is equal to y.
231
232
233        For small arguments 0 < y < exp(-2), the program computes
234        z = sqrt( -2.0 * log(y) );  then the approximation is
235        x = z - log(z)/z  - (1/z) P(1/z) / Q(1/z).
236        There are two rational functions P/Q, one for 0 < y < exp(-32)
237        and the other for y up to exp(-2).  For larger arguments,
238        w = y - 0.5, and  x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)).
239
240        ACCURACY:
241
242                             Relative error:
243        arithmetic   domain        # trials      peak         rms
244           IEEE     0.125, 1        20000       7.2e-16     1.3e-16
245           IEEE     3e-308, 0.135   50000       4.6e-16     9.8e-17
246
247        Cephes Math Library Release 2.8:  June, 2000
248        Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
249        *************************************************************************/
250        public static double invnormaldistribution(double y0)
251        {
252            double result = 0;
253            double expm2 = 0;
254            double s2pi = 0;
255            double x = 0;
256            double y = 0;
257            double z = 0;
258            double y2 = 0;
259            double x0 = 0;
260            double x1 = 0;
261            int code = 0;
262            double p0 = 0;
263            double q0 = 0;
264            double p1 = 0;
265            double q1 = 0;
266            double p2 = 0;
267            double q2 = 0;
268
269            expm2 = 0.13533528323661269189;
270            s2pi = 2.50662827463100050242;
271            if( y0<=0 )
272            {
273                result = -AP.Math.MaxRealNumber;
274                return result;
275            }
276            if( y0>=1 )
277            {
278                result = AP.Math.MaxRealNumber;
279                return result;
280            }
281            code = 1;
282            y = y0;
283            if( y>1.0-expm2 )
284            {
285                y = 1.0-y;
286                code = 0;
287            }
288            if( y>expm2 )
289            {
290                y = y-0.5;
291                y2 = y*y;
292                p0 = -59.9633501014107895267;
293                p0 = 98.0010754185999661536+y2*p0;
294                p0 = -56.6762857469070293439+y2*p0;
295                p0 = 13.9312609387279679503+y2*p0;
296                p0 = -1.23916583867381258016+y2*p0;
297                q0 = 1;
298                q0 = 1.95448858338141759834+y2*q0;
299                q0 = 4.67627912898881538453+y2*q0;
300                q0 = 86.3602421390890590575+y2*q0;
301                q0 = -225.462687854119370527+y2*q0;
302                q0 = 200.260212380060660359+y2*q0;
303                q0 = -82.0372256168333339912+y2*q0;
304                q0 = 15.9056225126211695515+y2*q0;
305                q0 = -1.18331621121330003142+y2*q0;
306                x = y+y*y2*p0/q0;
307                x = x*s2pi;
308                result = x;
309                return result;
310            }
311            x = Math.Sqrt(-(2.0*Math.Log(y)));
312            x0 = x-Math.Log(x)/x;
313            z = 1.0/x;
314            if( x<8.0 )
315            {
316                p1 = 4.05544892305962419923;
317                p1 = 31.5251094599893866154+z*p1;
318                p1 = 57.1628192246421288162+z*p1;
319                p1 = 44.0805073893200834700+z*p1;
320                p1 = 14.6849561928858024014+z*p1;
321                p1 = 2.18663306850790267539+z*p1;
322                p1 = -(1.40256079171354495875*0.1)+z*p1;
323                p1 = -(3.50424626827848203418*0.01)+z*p1;
324                p1 = -(8.57456785154685413611*0.0001)+z*p1;
325                q1 = 1;
326                q1 = 15.7799883256466749731+z*q1;
327                q1 = 45.3907635128879210584+z*q1;
328                q1 = 41.3172038254672030440+z*q1;
329                q1 = 15.0425385692907503408+z*q1;
330                q1 = 2.50464946208309415979+z*q1;
331                q1 = -(1.42182922854787788574*0.1)+z*q1;
332                q1 = -(3.80806407691578277194*0.01)+z*q1;
333                q1 = -(9.33259480895457427372*0.0001)+z*q1;
334                x1 = z*p1/q1;
335            }
336            else
337            {
338                p2 = 3.23774891776946035970;
339                p2 = 6.91522889068984211695+z*p2;
340                p2 = 3.93881025292474443415+z*p2;
341                p2 = 1.33303460815807542389+z*p2;
342                p2 = 2.01485389549179081538*0.1+z*p2;
343                p2 = 1.23716634817820021358*0.01+z*p2;
344                p2 = 3.01581553508235416007*0.0001+z*p2;
345                p2 = 2.65806974686737550832*0.000001+z*p2;
346                p2 = 6.23974539184983293730*0.000000001+z*p2;
347                q2 = 1;
348                q2 = 6.02427039364742014255+z*q2;
349                q2 = 3.67983563856160859403+z*q2;
350                q2 = 1.37702099489081330271+z*q2;
351                q2 = 2.16236993594496635890*0.1+z*q2;
352                q2 = 1.34204006088543189037*0.01+z*q2;
353                q2 = 3.28014464682127739104*0.0001+z*q2;
354                q2 = 2.89247864745380683936*0.000001+z*q2;
355                q2 = 6.79019408009981274425*0.000000001+z*q2;
356                x1 = z*p2/q2;
357            }
358            x = x0-x1;
359            if( code!=0 )
360            {
361                x = -x;
362            }
363            result = x;
364            return result;
365        }
366    }
367}
Note: See TracBrowser for help on using the repository browser.