ibetaf.cs

数理统计Stutent s检验源代码

                            {
                                di = di*di;
                            }
                            else
                            {
                                if( dir<-1 )
                                {
                                    di = 0.5*di;
                                }
                                else
                                {
                                    di = (yyy-y0)/(yh-yl);
                                }
                            }
                        }
                        dir = dir-1;
                    }
                    i = i+1;
                    mainlooppos = ihalvecycle;
                    continue;
                }
                else
                {
                    mainlooppos = breakihalvecycle;
                    continue;
                }
            }
            
            //
            // breakihalvecycle
            //
            if( mainlooppos==breakihalvecycle )
            {
                if( x0>=1.0 )
                {
                    x = 1.0-AP.Math.MachineEpsilon;
                    break;
                }
                if( x<=0.0 )
                {
                    x = 0.0;
                    break;
                }
                mainlooppos = newt;
                continue;
            }
            
            //
            // newt
            //
            if( mainlooppos==newt )
            {
                if( nflg!=0 )
                {
                    break;
                }
                nflg = 1;
                lgm = gammaf.lngamma(aaa+bbb, ref s)-gammaf.lngamma(aaa, ref s)-gammaf.lngamma(bbb, ref s);
                i = 0;
                mainlooppos = newtcycle;
                continue;
            }
            
            //
            // newtcycle
            //
            if( mainlooppos==newtcycle )
            {
                if( i<=7 )
                {
                    if( i!=0 )
                    {
                        yyy = incompletebeta(aaa, bbb, x);
                    }
                    if( yyy<yl )
                    {
                        x = x0;
                        yyy = yl;
                    }
                    else
                    {
                        if( yyy>yh )
                        {
                            x = x1;
                            yyy = yh;
                        }
                        else
                        {
                            if( yyy<y0 )
                            {
                                x0 = x;
                                yl = yyy;
                            }
                            else
                            {
                                x1 = x;
                                yh = yyy;
                            }
                        }
                    }
                    if( x==1.0 | x==0.0 )
                    {
                        mainlooppos = breaknewtcycle;
                        continue;
                    }
                    d = (aaa-1.0)*Math.Log(x)+(bbb-1.0)*Math.Log(1.0-x)+lgm;
                    if( d<Math.Log(AP.Math.MinRealNumber) )
                    {
                        break;
                    }
                    if( d>Math.Log(AP.Math.MaxRealNumber) )
                    {
                        mainlooppos = breaknewtcycle;
                        continue;
                    }
                    d = Math.Exp(d);
                    d = (yyy-y0)/d;
                    xt = x-d;
                    if( xt<=x0 )
                    {
                        yyy = (x-x0)/(x1-x0);
                        xt = x0+0.5*yyy*(x-x0);
                        if( xt<=0.0 )
                        {
                            mainlooppos = breaknewtcycle;
                            continue;
                        }
                    }
                    if( xt>=x1 )
                    {
                        yyy = (x1-x)/(x1-x0);
                        xt = x1-0.5*yyy*(x1-x);
                        if( xt>=1.0 )
                        {
                            mainlooppos = breaknewtcycle;
                            continue;
                        }
                    }
                    x = xt;
                    if( Math.Abs(d/x)<128.0*AP.Math.MachineEpsilon )
                    {
                        break;
                    }
                    i = i+1;
                    mainlooppos = newtcycle;
                    continue;
                }
                else
                {
                    mainlooppos = breaknewtcycle;
                    continue;
                }
            }
            
            //
            // breaknewtcycle
            //
            if( mainlooppos==breaknewtcycle )
            {
                dithresh = 256.0*AP.Math.MachineEpsilon;
                mainlooppos = ihalve;
                continue;
            }
        }
        
        //
        // done
        //
        if( rflg!=0 )
        {
            if( x<=AP.Math.MachineEpsilon )
            {
                x = 1.0-AP.Math.MachineEpsilon;
            }
            else
            {
                x = 1.0-x;
            }
        }
        result = x;
        return result;
    }


    /*************************************************************************
    Continued fraction expansion #1 for incomplete beta integral

    Cephes Math Library, Release 2.8:  June, 2000
    Copyright 1984, 1995, 2000 by Stephen L. Moshier
    *************************************************************************/
    private static double incompletebetafe(double a,
        double b,
        double x,
        double big,
        double biginv)
    {
        double result = 0;
        double xk = 0;
        double pk = 0;
        double pkm1 = 0;
        double pkm2 = 0;
        double qk = 0;
        double qkm1 = 0;
        double qkm2 = 0;
        double k1 = 0;
        double k2 = 0;
        double k3 = 0;
        double k4 = 0;
        double k5 = 0;
        double k6 = 0;
        double k7 = 0;
        double k8 = 0;
        double r = 0;
        double t = 0;
        double ans = 0;
        double thresh = 0;
        int n = 0;

        k1 = a;
        k2 = a+b;
        k3 = a;
        k4 = a+1.0;
        k5 = 1.0;
        k6 = b-1.0;
        k7 = k4;
        k8 = a+2.0;
        pkm2 = 0.0;
        qkm2 = 1.0;
        pkm1 = 1.0;
        qkm1 = 1.0;
        ans = 1.0;
        r = 1.0;
        n = 0;
        thresh = 3.0*AP.Math.MachineEpsilon;
        do
        {
            xk = -(x*k1*k2/(k3*k4));
            pk = pkm1+pkm2*xk;
            qk = qkm1+qkm2*xk;
            pkm2 = pkm1;
            pkm1 = pk;
            qkm2 = qkm1;
            qkm1 = qk;
            xk = x*k5*k6/(k7*k8);
            pk = pkm1+pkm2*xk;
            qk = qkm1+qkm2*xk;
            pkm2 = pkm1;
            pkm1 = pk;
            qkm2 = qkm1;
            qkm1 = qk;
            if( qk!=0 )
            {
                r = pk/qk;
            }
            if( r!=0 )
            {
                t = Math.Abs((ans-r)/r);
                ans = r;
            }
            else
            {
                t = 1.0;
            }
            if( t<thresh )
            {
                break;
            }
            k1 = k1+1.0;
            k2 = k2+1.0;
            k3 = k3+2.0;
            k4 = k4+2.0;
            k5 = k5+1.0;
            k6 = k6-1.0;
            k7 = k7+2.0;
            k8 = k8+2.0;
            if( Math.Abs(qk)+Math.Abs(pk)>big )
            {
                pkm2 = pkm2*biginv;
                pkm1 = pkm1*biginv;
                qkm2 = qkm2*biginv;
                qkm1 = qkm1*biginv;
            }
            if( Math.Abs(qk)<biginv | Math.Abs(pk)<biginv )
            {
                pkm2 = pkm2*big;
                pkm1 = pkm1*big;
                qkm2 = qkm2*big;
                qkm1 = qkm1*big;
            }
            n = n+1;
        }
        while( n!=300 );
        result = ans;
        return result;
    }


    /*************************************************************************
    Continued fraction expansion #2
    for incomplete beta integral

    Cephes Math Library, Release 2.8:  June, 2000
    Copyright 1984, 1995, 2000 by Stephen L. Moshier
    *************************************************************************/
    private static double incompletebetafe2(double a,
        double b,
        double x,
        double big,
        double biginv)
    {
        double result = 0;
        double xk = 0;
        double pk = 0;
        double pkm1 = 0;
        double pkm2 = 0;
        double qk = 0;
        double qkm1 = 0;
        double qkm2 = 0;
        double k1 = 0;
        double k2 = 0;
        double k3 = 0;
        double k4 = 0;
        double k5 = 0;
        double k6 = 0;
        double k7 = 0;
        double k8 = 0;
        double r = 0;
        double t = 0;
        double ans = 0;
        double z = 0;
        double thresh = 0;
        int n = 0;

        k1 = a;
        k2 = b-1.0;
        k3 = a;
        k4 = a+1.0;
        k5 = 1.0;
        k6 = a+b;
        k7 = a+1.0;
        k8 = a+2.0;
        pkm2 = 0.0;
        qkm2 = 1.0;
        pkm1 = 1.0;
        qkm1 = 1.0;
        z = x/(1.0-x);
        ans = 1.0;
        r = 1.0;
        n = 0;
        thresh = 3.0*AP.Math.MachineEpsilon;
        do
        {
            xk = -(z*k1*k2/(k3*k4));
            pk = pkm1+pkm2*xk;
            qk = qkm1+qkm2*xk;
            pkm2 = pkm1;
            pkm1 = pk;
            qkm2 = qkm1;
            qkm1 = qk;
            xk = z*k5*k6/(k7*k8);
            pk = pkm1+pkm2*xk;
            qk = qkm1+qkm2*xk;
            pkm2 = pkm1;
            pkm1 = pk;
            qkm2 = qkm1;
            qkm1 = qk;
            if( qk!=0 )
            {
                r = pk/qk;
            }
            if( r!=0 )
            {
                t = Math.Abs((ans-r)/r);
                ans = r;
            }
            else
            {
                t = 1.0;
            }
            if( t<thresh )
            {
                break;
            }
            k1 = k1+1.0;
            k2 = k2-1.0;
            k3 = k3+2.0;
            k4 = k4+2.0;
            k5 = k5+1.0;
            k6 = k6+1.0;
            k7 = k7+2.0;
            k8 = k8+2.0;
            if( Math.Abs(qk)+Math.Abs(pk)>big )
            {
                pkm2 = pkm2*biginv;
                pkm1 = pkm1*biginv;
                qkm2 = qkm2*biginv;
                qkm1 = qkm1*biginv;
            }
            if( Math.Abs(qk)<biginv | Math.Abs(pk)<biginv )
            {
                pkm2 = pkm2*big;
                pkm1 = pkm1*big;
                qkm2 = qkm2*big;
                qkm1 = qkm1*big;
            }
            n = n+1;
        }
        while( n!=300 );
        result = ans;
        return result;
    }


    /*************************************************************************
    Power series for incomplete beta integral.
    Use when b*x is small and x not too close to 1.

    Cephes Math Library, Release 2.8:  June, 2000
    Copyright 1984, 1995, 2000 by Stephen L. Moshier
    *************************************************************************/
    private static double incompletebetaps(double a,
        double b,
        double x,
        double maxgam)
    {
        double result = 0;
        double s = 0;
        double t = 0;
        double u = 0;
        double v = 0;
        double n = 0;
        double t1 = 0;
        double z = 0;
        double ai = 0;
        double sg = 0;

        ai = 1.0/a;
        u = (1.0-b)*x;
        v = u/(a+1.0);
        t1 = v;
        t = u;
        n = 2.0;
        s = 0.0;
        z = AP.Math.MachineEpsilon*ai;
        while( Math.Abs(v)>z )
        {
            u = (n-b)*x/n;
            t = t*u;
            v = t/(a+n);
            s = s+v;
            n = n+1.0;
        }
        s = s+t1;
        s = s+ai;
        u = a*Math.Log(x);
        if( a+b<maxgam & Math.Abs(u)<Math.Log(AP.Math.MaxRealNumber) )
        {
            t = gammaf.gamma(a+b)/(gammaf.gamma(a)*gammaf.gamma(b));
            s = s*t*Math.Pow(x, a);
        }
        else
        {
            t = gammaf.lngamma(a+b, ref sg)-gammaf.lngamma(a, ref sg)-gammaf.lngamma(b, ref sg)+u+Math.Log(s);
            if( t<Math.Log(AP.Math.MinRealNumber) )
            {
                s = 0.0;
            }
            else
            {
                s = Math.Exp(t);
            }
        }
        result = s;
        return result;
    }
}