Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/ALGLIB/ibetaf.cs @ 2575

Last change on this file since 2575 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: 32.0 KB
Line 
1/*************************************************************************
2Cephes Math Library Release 2.8:  June, 2000
3Copyright 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 ibetaf
33    {
34        /*************************************************************************
35        Incomplete beta integral
36
37        Returns incomplete beta integral of the arguments, evaluated
38        from zero to x.  The function is defined as
39
40                         x
41            -            -
42           | (a+b)      | |  a-1     b-1
43         -----------    |   t   (1-t)   dt.
44          -     -     | |
45         | (a) | (b)   -
46                        0
47
48        The domain of definition is 0 <= x <= 1.  In this
49        implementation a and b are restricted to positive values.
50        The integral from x to 1 may be obtained by the symmetry
51        relation
52
53           1 - incbet( a, b, x )  =  incbet( b, a, 1-x ).
54
55        The integral is evaluated by a continued fraction expansion
56        or, when b*x is small, by a power series.
57
58        ACCURACY:
59
60        Tested at uniformly distributed random points (a,b,x) with a and b
61        in "domain" and x between 0 and 1.
62                                               Relative error
63        arithmetic   domain     # trials      peak         rms
64           IEEE      0,5         10000       6.9e-15     4.5e-16
65           IEEE      0,85       250000       2.2e-13     1.7e-14
66           IEEE      0,1000      30000       5.3e-12     6.3e-13
67           IEEE      0,10000    250000       9.3e-11     7.1e-12
68           IEEE      0,100000    10000       8.7e-10     4.8e-11
69        Outputs smaller than the IEEE gradual underflow threshold
70        were excluded from these statistics.
71
72        Cephes Math Library, Release 2.8:  June, 2000
73        Copyright 1984, 1995, 2000 by Stephen L. Moshier
74        *************************************************************************/
75        public static double incompletebeta(double a,
76            double b,
77            double x)
78        {
79            double result = 0;
80            double t = 0;
81            double xc = 0;
82            double w = 0;
83            double y = 0;
84            int flag = 0;
85            double sg = 0;
86            double big = 0;
87            double biginv = 0;
88            double maxgam = 0;
89            double minlog = 0;
90            double maxlog = 0;
91
92            big = 4.503599627370496e15;
93            biginv = 2.22044604925031308085e-16;
94            maxgam = 171.624376956302725;
95            minlog = Math.Log(AP.Math.MinRealNumber);
96            maxlog = Math.Log(AP.Math.MaxRealNumber);
97            System.Diagnostics.Debug.Assert((double)(a)>(double)(0) & (double)(b)>(double)(0), "Domain error in IncompleteBeta");
98            System.Diagnostics.Debug.Assert((double)(x)>=(double)(0) & (double)(x)<=(double)(1), "Domain error in IncompleteBeta");
99            if( (double)(x)==(double)(0) )
100            {
101                result = 0;
102                return result;
103            }
104            if( (double)(x)==(double)(1) )
105            {
106                result = 1;
107                return result;
108            }
109            flag = 0;
110            if( (double)(b*x)<=(double)(1.0) & (double)(x)<=(double)(0.95) )
111            {
112                result = incompletebetaps(a, b, x, maxgam);
113                return result;
114            }
115            w = 1.0-x;
116            if( (double)(x)>(double)(a/(a+b)) )
117            {
118                flag = 1;
119                t = a;
120                a = b;
121                b = t;
122                xc = x;
123                x = w;
124            }
125            else
126            {
127                xc = w;
128            }
129            if( flag==1 & (double)(b*x)<=(double)(1.0) & (double)(x)<=(double)(0.95) )
130            {
131                t = incompletebetaps(a, b, x, maxgam);
132                if( (double)(t)<=(double)(AP.Math.MachineEpsilon) )
133                {
134                    result = 1.0-AP.Math.MachineEpsilon;
135                }
136                else
137                {
138                    result = 1.0-t;
139                }
140                return result;
141            }
142            y = x*(a+b-2.0)-(a-1.0);
143            if( (double)(y)<(double)(0.0) )
144            {
145                w = incompletebetafe(a, b, x, big, biginv);
146            }
147            else
148            {
149                w = incompletebetafe2(a, b, x, big, biginv)/xc;
150            }
151            y = a*Math.Log(x);
152            t = b*Math.Log(xc);
153            if( (double)(a+b)<(double)(maxgam) & (double)(Math.Abs(y))<(double)(maxlog) & (double)(Math.Abs(t))<(double)(maxlog) )
154            {
155                t = Math.Pow(xc, b);
156                t = t*Math.Pow(x, a);
157                t = t/a;
158                t = t*w;
159                t = t*(gammafunc.gamma(a+b)/(gammafunc.gamma(a)*gammafunc.gamma(b)));
160                if( flag==1 )
161                {
162                    if( (double)(t)<=(double)(AP.Math.MachineEpsilon) )
163                    {
164                        result = 1.0-AP.Math.MachineEpsilon;
165                    }
166                    else
167                    {
168                        result = 1.0-t;
169                    }
170                }
171                else
172                {
173                    result = t;
174                }
175                return result;
176            }
177            y = y+t+gammafunc.lngamma(a+b, ref sg)-gammafunc.lngamma(a, ref sg)-gammafunc.lngamma(b, ref sg);
178            y = y+Math.Log(w/a);
179            if( (double)(y)<(double)(minlog) )
180            {
181                t = 0.0;
182            }
183            else
184            {
185                t = Math.Exp(y);
186            }
187            if( flag==1 )
188            {
189                if( (double)(t)<=(double)(AP.Math.MachineEpsilon) )
190                {
191                    t = 1.0-AP.Math.MachineEpsilon;
192                }
193                else
194                {
195                    t = 1.0-t;
196                }
197            }
198            result = t;
199            return result;
200        }
201
202
203        /*************************************************************************
204        Inverse of imcomplete beta integral
205
206        Given y, the function finds x such that
207
208         incbet( a, b, x ) = y .
209
210        The routine performs interval halving or Newton iterations to find the
211        root of incbet(a,b,x) - y = 0.
212
213
214        ACCURACY:
215
216                             Relative error:
217                       x     a,b
218        arithmetic   domain  domain  # trials    peak       rms
219           IEEE      0,1    .5,10000   50000    5.8e-12   1.3e-13
220           IEEE      0,1   .25,100    100000    1.8e-13   3.9e-15
221           IEEE      0,1     0,5       50000    1.1e-12   5.5e-15
222        With a and b constrained to half-integer or integer values:
223           IEEE      0,1    .5,10000   50000    5.8e-12   1.1e-13
224           IEEE      0,1    .5,100    100000    1.7e-14   7.9e-16
225        With a = .5, b constrained to half-integer or integer values:
226           IEEE      0,1    .5,10000   10000    8.3e-11   1.0e-11
227
228        Cephes Math Library Release 2.8:  June, 2000
229        Copyright 1984, 1996, 2000 by Stephen L. Moshier
230        *************************************************************************/
231        public static double invincompletebeta(double a,
232            double b,
233            double y)
234        {
235            double result = 0;
236            double aaa = 0;
237            double bbb = 0;
238            double y0 = 0;
239            double d = 0;
240            double yyy = 0;
241            double x = 0;
242            double x0 = 0;
243            double x1 = 0;
244            double lgm = 0;
245            double yp = 0;
246            double di = 0;
247            double dithresh = 0;
248            double yl = 0;
249            double yh = 0;
250            double xt = 0;
251            int i = 0;
252            int rflg = 0;
253            int dir = 0;
254            int nflg = 0;
255            double s = 0;
256            int mainlooppos = 0;
257            int ihalve = 0;
258            int ihalvecycle = 0;
259            int newt = 0;
260            int newtcycle = 0;
261            int breaknewtcycle = 0;
262            int breakihalvecycle = 0;
263
264            i = 0;
265            System.Diagnostics.Debug.Assert((double)(y)>=(double)(0) & (double)(y)<=(double)(1), "Domain error in InvIncompleteBeta");
266            if( (double)(y)==(double)(0) )
267            {
268                result = 0;
269                return result;
270            }
271            if( (double)(y)==(double)(1.0) )
272            {
273                result = 1;
274                return result;
275            }
276            x0 = 0.0;
277            yl = 0.0;
278            x1 = 1.0;
279            yh = 1.0;
280            nflg = 0;
281            mainlooppos = 0;
282            ihalve = 1;
283            ihalvecycle = 2;
284            newt = 3;
285            newtcycle = 4;
286            breaknewtcycle = 5;
287            breakihalvecycle = 6;
288            while( true )
289            {
290               
291                //
292                // start
293                //
294                if( mainlooppos==0 )
295                {
296                    if( (double)(a)<=(double)(1.0) | (double)(b)<=(double)(1.0) )
297                    {
298                        dithresh = 1.0e-6;
299                        rflg = 0;
300                        aaa = a;
301                        bbb = b;
302                        y0 = y;
303                        x = aaa/(aaa+bbb);
304                        yyy = incompletebeta(aaa, bbb, x);
305                        mainlooppos = ihalve;
306                        continue;
307                    }
308                    else
309                    {
310                        dithresh = 1.0e-4;
311                    }
312                    yp = -normaldistr.invnormaldistribution(y);
313                    if( (double)(y)>(double)(0.5) )
314                    {
315                        rflg = 1;
316                        aaa = b;
317                        bbb = a;
318                        y0 = 1.0-y;
319                        yp = -yp;
320                    }
321                    else
322                    {
323                        rflg = 0;
324                        aaa = a;
325                        bbb = b;
326                        y0 = y;
327                    }
328                    lgm = (yp*yp-3.0)/6.0;
329                    x = 2.0/(1.0/(2.0*aaa-1.0)+1.0/(2.0*bbb-1.0));
330                    d = yp*Math.Sqrt(x+lgm)/x-(1.0/(2.0*bbb-1.0)-1.0/(2.0*aaa-1.0))*(lgm+5.0/6.0-2.0/(3.0*x));
331                    d = 2.0*d;
332                    if( (double)(d)<(double)(Math.Log(AP.Math.MinRealNumber)) )
333                    {
334                        x = 0;
335                        break;
336                    }
337                    x = aaa/(aaa+bbb*Math.Exp(d));
338                    yyy = incompletebeta(aaa, bbb, x);
339                    yp = (yyy-y0)/y0;
340                    if( (double)(Math.Abs(yp))<(double)(0.2) )
341                    {
342                        mainlooppos = newt;
343                        continue;
344                    }
345                    mainlooppos = ihalve;
346                    continue;
347                }
348               
349                //
350                // ihalve
351                //
352                if( mainlooppos==ihalve )
353                {
354                    dir = 0;
355                    di = 0.5;
356                    i = 0;
357                    mainlooppos = ihalvecycle;
358                    continue;
359                }
360               
361                //
362                // ihalvecycle
363                //
364                if( mainlooppos==ihalvecycle )
365                {
366                    if( i<=99 )
367                    {
368                        if( i!=0 )
369                        {
370                            x = x0+di*(x1-x0);
371                            if( (double)(x)==(double)(1.0) )
372                            {
373                                x = 1.0-AP.Math.MachineEpsilon;
374                            }
375                            if( (double)(x)==(double)(0.0) )
376                            {
377                                di = 0.5;
378                                x = x0+di*(x1-x0);
379                                if( (double)(x)==(double)(0.0) )
380                                {
381                                    break;
382                                }
383                            }
384                            yyy = incompletebeta(aaa, bbb, x);
385                            yp = (x1-x0)/(x1+x0);
386                            if( (double)(Math.Abs(yp))<(double)(dithresh) )
387                            {
388                                mainlooppos = newt;
389                                continue;
390                            }
391                            yp = (yyy-y0)/y0;
392                            if( (double)(Math.Abs(yp))<(double)(dithresh) )
393                            {
394                                mainlooppos = newt;
395                                continue;
396                            }
397                        }
398                        if( (double)(yyy)<(double)(y0) )
399                        {
400                            x0 = x;
401                            yl = yyy;
402                            if( dir<0 )
403                            {
404                                dir = 0;
405                                di = 0.5;
406                            }
407                            else
408                            {
409                                if( dir>3 )
410                                {
411                                    di = 1.0-(1.0-di)*(1.0-di);
412                                }
413                                else
414                                {
415                                    if( dir>1 )
416                                    {
417                                        di = 0.5*di+0.5;
418                                    }
419                                    else
420                                    {
421                                        di = (y0-yyy)/(yh-yl);
422                                    }
423                                }
424                            }
425                            dir = dir+1;
426                            if( (double)(x0)>(double)(0.75) )
427                            {
428                                if( rflg==1 )
429                                {
430                                    rflg = 0;
431                                    aaa = a;
432                                    bbb = b;
433                                    y0 = y;
434                                }
435                                else
436                                {
437                                    rflg = 1;
438                                    aaa = b;
439                                    bbb = a;
440                                    y0 = 1.0-y;
441                                }
442                                x = 1.0-x;
443                                yyy = incompletebeta(aaa, bbb, x);
444                                x0 = 0.0;
445                                yl = 0.0;
446                                x1 = 1.0;
447                                yh = 1.0;
448                                mainlooppos = ihalve;
449                                continue;
450                            }
451                        }
452                        else
453                        {
454                            x1 = x;
455                            if( rflg==1 & (double)(x1)<(double)(AP.Math.MachineEpsilon) )
456                            {
457                                x = 0.0;
458                                break;
459                            }
460                            yh = yyy;
461                            if( dir>0 )
462                            {
463                                dir = 0;
464                                di = 0.5;
465                            }
466                            else
467                            {
468                                if( dir<-3 )
469                                {
470                                    di = di*di;
471                                }
472                                else
473                                {
474                                    if( dir<-1 )
475                                    {
476                                        di = 0.5*di;
477                                    }
478                                    else
479                                    {
480                                        di = (yyy-y0)/(yh-yl);
481                                    }
482                                }
483                            }
484                            dir = dir-1;
485                        }
486                        i = i+1;
487                        mainlooppos = ihalvecycle;
488                        continue;
489                    }
490                    else
491                    {
492                        mainlooppos = breakihalvecycle;
493                        continue;
494                    }
495                }
496               
497                //
498                // breakihalvecycle
499                //
500                if( mainlooppos==breakihalvecycle )
501                {
502                    if( (double)(x0)>=(double)(1.0) )
503                    {
504                        x = 1.0-AP.Math.MachineEpsilon;
505                        break;
506                    }
507                    if( (double)(x)<=(double)(0.0) )
508                    {
509                        x = 0.0;
510                        break;
511                    }
512                    mainlooppos = newt;
513                    continue;
514                }
515               
516                //
517                // newt
518                //
519                if( mainlooppos==newt )
520                {
521                    if( nflg!=0 )
522                    {
523                        break;
524                    }
525                    nflg = 1;
526                    lgm = gammafunc.lngamma(aaa+bbb, ref s)-gammafunc.lngamma(aaa, ref s)-gammafunc.lngamma(bbb, ref s);
527                    i = 0;
528                    mainlooppos = newtcycle;
529                    continue;
530                }
531               
532                //
533                // newtcycle
534                //
535                if( mainlooppos==newtcycle )
536                {
537                    if( i<=7 )
538                    {
539                        if( i!=0 )
540                        {
541                            yyy = incompletebeta(aaa, bbb, x);
542                        }
543                        if( (double)(yyy)<(double)(yl) )
544                        {
545                            x = x0;
546                            yyy = yl;
547                        }
548                        else
549                        {
550                            if( (double)(yyy)>(double)(yh) )
551                            {
552                                x = x1;
553                                yyy = yh;
554                            }
555                            else
556                            {
557                                if( (double)(yyy)<(double)(y0) )
558                                {
559                                    x0 = x;
560                                    yl = yyy;
561                                }
562                                else
563                                {
564                                    x1 = x;
565                                    yh = yyy;
566                                }
567                            }
568                        }
569                        if( (double)(x)==(double)(1.0) | (double)(x)==(double)(0.0) )
570                        {
571                            mainlooppos = breaknewtcycle;
572                            continue;
573                        }
574                        d = (aaa-1.0)*Math.Log(x)+(bbb-1.0)*Math.Log(1.0-x)+lgm;
575                        if( (double)(d)<(double)(Math.Log(AP.Math.MinRealNumber)) )
576                        {
577                            break;
578                        }
579                        if( (double)(d)>(double)(Math.Log(AP.Math.MaxRealNumber)) )
580                        {
581                            mainlooppos = breaknewtcycle;
582                            continue;
583                        }
584                        d = Math.Exp(d);
585                        d = (yyy-y0)/d;
586                        xt = x-d;
587                        if( (double)(xt)<=(double)(x0) )
588                        {
589                            yyy = (x-x0)/(x1-x0);
590                            xt = x0+0.5*yyy*(x-x0);
591                            if( (double)(xt)<=(double)(0.0) )
592                            {
593                                mainlooppos = breaknewtcycle;
594                                continue;
595                            }
596                        }
597                        if( (double)(xt)>=(double)(x1) )
598                        {
599                            yyy = (x1-x)/(x1-x0);
600                            xt = x1-0.5*yyy*(x1-x);
601                            if( (double)(xt)>=(double)(1.0) )
602                            {
603                                mainlooppos = breaknewtcycle;
604                                continue;
605                            }
606                        }
607                        x = xt;
608                        if( (double)(Math.Abs(d/x))<(double)(128.0*AP.Math.MachineEpsilon) )
609                        {
610                            break;
611                        }
612                        i = i+1;
613                        mainlooppos = newtcycle;
614                        continue;
615                    }
616                    else
617                    {
618                        mainlooppos = breaknewtcycle;
619                        continue;
620                    }
621                }
622               
623                //
624                // breaknewtcycle
625                //
626                if( mainlooppos==breaknewtcycle )
627                {
628                    dithresh = 256.0*AP.Math.MachineEpsilon;
629                    mainlooppos = ihalve;
630                    continue;
631                }
632            }
633           
634            //
635            // done
636            //
637            if( rflg!=0 )
638            {
639                if( (double)(x)<=(double)(AP.Math.MachineEpsilon) )
640                {
641                    x = 1.0-AP.Math.MachineEpsilon;
642                }
643                else
644                {
645                    x = 1.0-x;
646                }
647            }
648            result = x;
649            return result;
650        }
651
652
653        /*************************************************************************
654        Continued fraction expansion #1 for incomplete beta integral
655
656        Cephes Math Library, Release 2.8:  June, 2000
657        Copyright 1984, 1995, 2000 by Stephen L. Moshier
658        *************************************************************************/
659        private static double incompletebetafe(double a,
660            double b,
661            double x,
662            double big,
663            double biginv)
664        {
665            double result = 0;
666            double xk = 0;
667            double pk = 0;
668            double pkm1 = 0;
669            double pkm2 = 0;
670            double qk = 0;
671            double qkm1 = 0;
672            double qkm2 = 0;
673            double k1 = 0;
674            double k2 = 0;
675            double k3 = 0;
676            double k4 = 0;
677            double k5 = 0;
678            double k6 = 0;
679            double k7 = 0;
680            double k8 = 0;
681            double r = 0;
682            double t = 0;
683            double ans = 0;
684            double thresh = 0;
685            int n = 0;
686
687            k1 = a;
688            k2 = a+b;
689            k3 = a;
690            k4 = a+1.0;
691            k5 = 1.0;
692            k6 = b-1.0;
693            k7 = k4;
694            k8 = a+2.0;
695            pkm2 = 0.0;
696            qkm2 = 1.0;
697            pkm1 = 1.0;
698            qkm1 = 1.0;
699            ans = 1.0;
700            r = 1.0;
701            n = 0;
702            thresh = 3.0*AP.Math.MachineEpsilon;
703            do
704            {
705                xk = -(x*k1*k2/(k3*k4));
706                pk = pkm1+pkm2*xk;
707                qk = qkm1+qkm2*xk;
708                pkm2 = pkm1;
709                pkm1 = pk;
710                qkm2 = qkm1;
711                qkm1 = qk;
712                xk = x*k5*k6/(k7*k8);
713                pk = pkm1+pkm2*xk;
714                qk = qkm1+qkm2*xk;
715                pkm2 = pkm1;
716                pkm1 = pk;
717                qkm2 = qkm1;
718                qkm1 = qk;
719                if( (double)(qk)!=(double)(0) )
720                {
721                    r = pk/qk;
722                }
723                if( (double)(r)!=(double)(0) )
724                {
725                    t = Math.Abs((ans-r)/r);
726                    ans = r;
727                }
728                else
729                {
730                    t = 1.0;
731                }
732                if( (double)(t)<(double)(thresh) )
733                {
734                    break;
735                }
736                k1 = k1+1.0;
737                k2 = k2+1.0;
738                k3 = k3+2.0;
739                k4 = k4+2.0;
740                k5 = k5+1.0;
741                k6 = k6-1.0;
742                k7 = k7+2.0;
743                k8 = k8+2.0;
744                if( (double)(Math.Abs(qk)+Math.Abs(pk))>(double)(big) )
745                {
746                    pkm2 = pkm2*biginv;
747                    pkm1 = pkm1*biginv;
748                    qkm2 = qkm2*biginv;
749                    qkm1 = qkm1*biginv;
750                }
751                if( (double)(Math.Abs(qk))<(double)(biginv) | (double)(Math.Abs(pk))<(double)(biginv) )
752                {
753                    pkm2 = pkm2*big;
754                    pkm1 = pkm1*big;
755                    qkm2 = qkm2*big;
756                    qkm1 = qkm1*big;
757                }
758                n = n+1;
759            }
760            while( n!=300 );
761            result = ans;
762            return result;
763        }
764
765
766        /*************************************************************************
767        Continued fraction expansion #2
768        for incomplete beta integral
769
770        Cephes Math Library, Release 2.8:  June, 2000
771        Copyright 1984, 1995, 2000 by Stephen L. Moshier
772        *************************************************************************/
773        private static double incompletebetafe2(double a,
774            double b,
775            double x,
776            double big,
777            double biginv)
778        {
779            double result = 0;
780            double xk = 0;
781            double pk = 0;
782            double pkm1 = 0;
783            double pkm2 = 0;
784            double qk = 0;
785            double qkm1 = 0;
786            double qkm2 = 0;
787            double k1 = 0;
788            double k2 = 0;
789            double k3 = 0;
790            double k4 = 0;
791            double k5 = 0;
792            double k6 = 0;
793            double k7 = 0;
794            double k8 = 0;
795            double r = 0;
796            double t = 0;
797            double ans = 0;
798            double z = 0;
799            double thresh = 0;
800            int n = 0;
801
802            k1 = a;
803            k2 = b-1.0;
804            k3 = a;
805            k4 = a+1.0;
806            k5 = 1.0;
807            k6 = a+b;
808            k7 = a+1.0;
809            k8 = a+2.0;
810            pkm2 = 0.0;
811            qkm2 = 1.0;
812            pkm1 = 1.0;
813            qkm1 = 1.0;
814            z = x/(1.0-x);
815            ans = 1.0;
816            r = 1.0;
817            n = 0;
818            thresh = 3.0*AP.Math.MachineEpsilon;
819            do
820            {
821                xk = -(z*k1*k2/(k3*k4));
822                pk = pkm1+pkm2*xk;
823                qk = qkm1+qkm2*xk;
824                pkm2 = pkm1;
825                pkm1 = pk;
826                qkm2 = qkm1;
827                qkm1 = qk;
828                xk = z*k5*k6/(k7*k8);
829                pk = pkm1+pkm2*xk;
830                qk = qkm1+qkm2*xk;
831                pkm2 = pkm1;
832                pkm1 = pk;
833                qkm2 = qkm1;
834                qkm1 = qk;
835                if( (double)(qk)!=(double)(0) )
836                {
837                    r = pk/qk;
838                }
839                if( (double)(r)!=(double)(0) )
840                {
841                    t = Math.Abs((ans-r)/r);
842                    ans = r;
843                }
844                else
845                {
846                    t = 1.0;
847                }
848                if( (double)(t)<(double)(thresh) )
849                {
850                    break;
851                }
852                k1 = k1+1.0;
853                k2 = k2-1.0;
854                k3 = k3+2.0;
855                k4 = k4+2.0;
856                k5 = k5+1.0;
857                k6 = k6+1.0;
858                k7 = k7+2.0;
859                k8 = k8+2.0;
860                if( (double)(Math.Abs(qk)+Math.Abs(pk))>(double)(big) )
861                {
862                    pkm2 = pkm2*biginv;
863                    pkm1 = pkm1*biginv;
864                    qkm2 = qkm2*biginv;
865                    qkm1 = qkm1*biginv;
866                }
867                if( (double)(Math.Abs(qk))<(double)(biginv) | (double)(Math.Abs(pk))<(double)(biginv) )
868                {
869                    pkm2 = pkm2*big;
870                    pkm1 = pkm1*big;
871                    qkm2 = qkm2*big;
872                    qkm1 = qkm1*big;
873                }
874                n = n+1;
875            }
876            while( n!=300 );
877            result = ans;
878            return result;
879        }
880
881
882        /*************************************************************************
883        Power series for incomplete beta integral.
884        Use when b*x is small and x not too close to 1.
885
886        Cephes Math Library, Release 2.8:  June, 2000
887        Copyright 1984, 1995, 2000 by Stephen L. Moshier
888        *************************************************************************/
889        private static double incompletebetaps(double a,
890            double b,
891            double x,
892            double maxgam)
893        {
894            double result = 0;
895            double s = 0;
896            double t = 0;
897            double u = 0;
898            double v = 0;
899            double n = 0;
900            double t1 = 0;
901            double z = 0;
902            double ai = 0;
903            double sg = 0;
904
905            ai = 1.0/a;
906            u = (1.0-b)*x;
907            v = u/(a+1.0);
908            t1 = v;
909            t = u;
910            n = 2.0;
911            s = 0.0;
912            z = AP.Math.MachineEpsilon*ai;
913            while( (double)(Math.Abs(v))>(double)(z) )
914            {
915                u = (n-b)*x/n;
916                t = t*u;
917                v = t/(a+n);
918                s = s+v;
919                n = n+1.0;
920            }
921            s = s+t1;
922            s = s+ai;
923            u = a*Math.Log(x);
924            if( (double)(a+b)<(double)(maxgam) & (double)(Math.Abs(u))<(double)(Math.Log(AP.Math.MaxRealNumber)) )
925            {
926                t = gammafunc.gamma(a+b)/(gammafunc.gamma(a)*gammafunc.gamma(b));
927                s = s*t*Math.Pow(x, a);
928            }
929            else
930            {
931                t = gammafunc.lngamma(a+b, ref sg)-gammafunc.lngamma(a, ref sg)-gammafunc.lngamma(b, ref sg)+u+Math.Log(s);
932                if( (double)(t)<(double)(Math.Log(AP.Math.MinRealNumber)) )
933                {
934                    s = 0.0;
935                }
936                else
937                {
938                    s = Math.Exp(t);
939                }
940            }
941            result = s;
942            return result;
943        }
944    }
945}
Note: See TracBrowser for help on using the repository browser.