dcdflib.c

unix环境下计算符合常用概率分布的随机变量的集合

    goto S10;
S20:
/*
                 TERMINATION
*/
    bfrac *= r;
    return bfrac;
}
void bgrat(double *a,double *b,double *x,double *y,double *w,
	   double *eps,int *ierr)
/*
-----------------------------------------------------------------------
     ASYMPTOTIC EXPANSION FOR IX(A,B) WHEN A IS LARGER THAN B.
     THE RESULT OF THE EXPANSION IS ADDED TO W. IT IS ASSUMED
     THAT A .GE. 15 AND B .LE. 1.  EPS IS THE TOLERANCE USED.
     IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS.
-----------------------------------------------------------------------
*/
{
static double bm1,bp2n,cn,coef,dj,j,l,lnx,n2,nu,p,q,r,s,sum,t,t2,u,v,z;
static int i,n,nm1;
static double c[30],d[30],T1;
/*
     ..
     .. Executable Statements ..
*/
    bm1 = *b-0.5e0-0.5e0;
    nu = *a+0.5e0*bm1;
    if(*y > 0.375e0) goto S10;
    T1 = -*y;
    lnx = alnrel(&T1);
    goto S20;
S10:
    lnx = log(*x);
S20:
    z = -(nu*lnx);
    if(*b*z == 0.0e0) goto S70;
/*
                 COMPUTATION OF THE EXPANSION
                 SET R = EXP(-Z)*Z**B/GAMMA(B)
*/
    r = *b*(1.0e0+gam1(b))*exp(*b*log(z));
    r *= (exp(*a*lnx)*exp(0.5e0*bm1*lnx));
    u = algdiv(b,a)+*b*log(nu);
    u = r*exp(-u);
    if(u == 0.0e0) goto S70;
    grat1(b,&z,&r,&p,&q,eps);
    v = 0.25e0*pow(1.0e0/nu,2.0);
    t2 = 0.25e0*lnx*lnx;
    l = *w/u;
    j = q/r;
    sum = j;
    t = cn = 1.0e0;
    n2 = 0.0e0;
    for(n=1; n<=30; n++) {
        bp2n = *b+n2;
        j = (bp2n*(bp2n+1.0e0)*j+(z+bp2n+1.0e0)*t)*v;
        n2 += 2.0e0;
        t *= t2;
        cn /= (n2*(n2+1.0e0));
        c[n-1] = cn;
        s = 0.0e0;
        if(n == 1) goto S40;
        nm1 = n-1;
        coef = *b-(double)n;
        for(i=1; i<=nm1; i++) {
            s += (coef*c[i-1]*d[n-i-1]);
            coef += *b;
        }
S40:
        d[n-1] = bm1*cn+s/(double)n;
        dj = d[n-1]*j;
        sum += dj;
        if(sum <= 0.0e0) goto S70;
        if(fabs(dj) <= *eps*(sum+l)) goto S60;
    }
S60:
/*
                    ADD THE RESULTS TO W
*/
    *ierr = 0;
    *w += (u*sum);
    return;
S70:
/*
               THE EXPANSION CANNOT BE COMPUTED
*/
    *ierr = 1;
    return;
}
double bpser(double *a,double *b,double *x,double *eps)
/*
-----------------------------------------------------------------------
     POWER SERIES EXPANSION FOR EVALUATING IX(A,B) WHEN B .LE. 1
     OR B*X .LE. 0.7.  EPS IS THE TOLERANCE USED.
-----------------------------------------------------------------------
*/
{
static double bpser,a0,apb,b0,c,n,sum,t,tol,u,w,z;
static int i,m;
/*
     ..
     .. Executable Statements ..
*/
    bpser = 0.0e0;
    if(*x == 0.0e0) return bpser;
/*
-----------------------------------------------------------------------
            COMPUTE THE FACTOR X**A/(A*BETA(A,B))
-----------------------------------------------------------------------
*/
    a0 = fifdmin1(*a,*b);
    if(a0 < 1.0e0) goto S10;
    z = *a*log(*x)-betaln(a,b);
    bpser = exp(z)/ *a;
    goto S100;
S10:
    b0 = fifdmax1(*a,*b);
    if(b0 >= 8.0e0) goto S90;
    if(b0 > 1.0e0) goto S40;
/*
            PROCEDURE FOR A0 .LT. 1 AND B0 .LE. 1
*/
    bpser = pow(*x,*a);
    if(bpser == 0.0e0) return bpser;
    apb = *a+*b;
    if(apb > 1.0e0) goto S20;
    z = 1.0e0+gam1(&apb);
    goto S30;
S20:
    u = *a+*b-1.e0;
    z = (1.0e0+gam1(&u))/apb;
S30:
    c = (1.0e0+gam1(a))*(1.0e0+gam1(b))/z;
    bpser *= (c*(*b/apb));
    goto S100;
S40:
/*
         PROCEDURE FOR A0 .LT. 1 AND 1 .LT. B0 .LT. 8
*/
    u = gamln1(&a0);
    m = b0-1.0e0;
    if(m < 1) goto S60;
    c = 1.0e0;
    for(i=1; i<=m; i++) {
        b0 -= 1.0e0;
        c *= (b0/(a0+b0));
    }
    u = log(c)+u;
S60:
    z = *a*log(*x)-u;
    b0 -= 1.0e0;
    apb = a0+b0;
    if(apb > 1.0e0) goto S70;
    t = 1.0e0+gam1(&apb);
    goto S80;
S70:
    u = a0+b0-1.e0;
    t = (1.0e0+gam1(&u))/apb;
S80:
    bpser = exp(z)*(a0/ *a)*(1.0e0+gam1(&b0))/t;
    goto S100;
S90:
/*
            PROCEDURE FOR A0 .LT. 1 AND B0 .GE. 8
*/
    u = gamln1(&a0)+algdiv(&a0,&b0);
    z = *a*log(*x)-u;
    bpser = a0/ *a*exp(z);
S100:
    if(bpser == 0.0e0 || *a <= 0.1e0**eps) return bpser;
/*
-----------------------------------------------------------------------
                     COMPUTE THE SERIES
-----------------------------------------------------------------------
*/
    sum = n = 0.0e0;
    c = 1.0e0;
    tol = *eps/ *a;
S110:
    n += 1.0e0;
    c *= ((0.5e0+(0.5e0-*b/n))**x);
    w = c/(*a+n);
    sum += w;
    if(fabs(w) > tol) goto S110;
    bpser *= (1.0e0+*a*sum);
    return bpser;
}
void bratio(double *a,double *b,double *x,double *y,double *w,
	    double *w1,int *ierr)
/*
-----------------------------------------------------------------------
 
            EVALUATION OF THE INCOMPLETE BETA FUNCTION IX(A,B)
 
                     --------------------
 
     IT IS ASSUMED THAT A AND B ARE NONNEGATIVE, AND THAT X .LE. 1
     AND Y = 1 - X.  BRATIO ASSIGNS W AND W1 THE VALUES
 
                      W  = IX(A,B)
                      W1 = 1 - IX(A,B)
 
     IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS.
     IF NO INPUT ERRORS ARE DETECTED THEN IERR IS SET TO 0 AND
     W AND W1 ARE COMPUTED. OTHERWISE, IF AN ERROR IS DETECTED,
     THEN W AND W1 ARE ASSIGNED THE VALUE 0 AND IERR IS SET TO
     ONE OF THE FOLLOWING VALUES ...
 
        IERR = 1  IF A OR B IS NEGATIVE
        IERR = 2  IF A = B = 0
        IERR = 3  IF X .LT. 0 OR X .GT. 1
        IERR = 4  IF Y .LT. 0 OR Y .GT. 1
        IERR = 5  IF X + Y .NE. 1
        IERR = 6  IF X = A = 0
        IERR = 7  IF Y = B = 0
 
--------------------
     WRITTEN BY ALFRED H. MORRIS, JR.
        NAVAL SURFACE WARFARE CENTER
        DAHLGREN, VIRGINIA
     REVISED ... NOV 1991
-----------------------------------------------------------------------
*/
{
static int K1 = 1;
static double a0,b0,eps,lambda,t,x0,y0,z;
static int ierr1,ind,n;
static double T2,T3,T4,T5;
/*
     ..
     .. Executable Statements ..
*/
/*
     ****** EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE SMALLEST
            FLOATING POINT NUMBER FOR WHICH 1.0 + EPS .GT. 1.0
*/
    eps = spmpar(&K1);
    *w = *w1 = 0.0e0;
    if(*a < 0.0e0 || *b < 0.0e0) goto S270;
    if(*a == 0.0e0 && *b == 0.0e0) goto S280;
    if(*x < 0.0e0 || *x > 1.0e0) goto S290;
    if(*y < 0.0e0 || *y > 1.0e0) goto S300;
    z = *x+*y-0.5e0-0.5e0;
    if(fabs(z) > 3.0e0*eps) goto S310;
    *ierr = 0;
    if(*x == 0.0e0) goto S210;
    if(*y == 0.0e0) goto S230;
    if(*a == 0.0e0) goto S240;
    if(*b == 0.0e0) goto S220;
    eps = fifdmax1(eps,1.e-15);
    if(fifdmax1(*a,*b) < 1.e-3*eps) goto S260;
    ind = 0;
    a0 = *a;
    b0 = *b;
    x0 = *x;
    y0 = *y;
    if(fifdmin1(a0,b0) > 1.0e0) goto S40;
/*
             PROCEDURE FOR A0 .LE. 1 OR B0 .LE. 1
*/
    if(*x <= 0.5e0) goto S10;
    ind = 1;
    a0 = *b;
    b0 = *a;
    x0 = *y;
    y0 = *x;
S10:
    if(b0 < fifdmin1(eps,eps*a0)) goto S90;
    if(a0 < fifdmin1(eps,eps*b0) && b0*x0 <= 1.0e0) goto S100;
    if(fifdmax1(a0,b0) > 1.0e0) goto S20;
    if(a0 >= fifdmin1(0.2e0,b0)) goto S110;
    if(pow(x0,a0) <= 0.9e0) goto S110;
    if(x0 >= 0.3e0) goto S120;
    n = 20;
    goto S140;
S20:
    if(b0 <= 1.0e0) goto S110;
    if(x0 >= 0.3e0) goto S120;
    if(x0 >= 0.1e0) goto S30;
    if(pow(x0*b0,a0) <= 0.7e0) goto S110;
S30:
    if(b0 > 15.0e0) goto S150;
    n = 20;
    goto S140;
S40:
/*
             PROCEDURE FOR A0 .GT. 1 AND B0 .GT. 1
*/
    if(*a > *b) goto S50;
    lambda = *a-(*a+*b)**x;
    goto S60;
S50:
    lambda = (*a+*b)**y-*b;
S60:
    if(lambda >= 0.0e0) goto S70;
    ind = 1;
    a0 = *b;
    b0 = *a;
    x0 = *y;
    y0 = *x;
    lambda = fabs(lambda);
S70:
    if(b0 < 40.0e0 && b0*x0 <= 0.7e0) goto S110;
    if(b0 < 40.0e0) goto S160;
    if(a0 > b0) goto S80;
    if(a0 <= 100.0e0) goto S130;
    if(lambda > 0.03e0*a0) goto S130;
    goto S200;
S80:
    if(b0 <= 100.0e0) goto S130;
    if(lambda > 0.03e0*b0) goto S130;
    goto S200;
S90:
/*
            EVALUATION OF THE APPROPRIATE ALGORITHM
*/
    *w = fpser(&a0,&b0,&x0,&eps);
    *w1 = 0.5e0+(0.5e0-*w);
    goto S250;
S100:
    *w1 = apser(&a0,&b0,&x0,&eps);
    *w = 0.5e0+(0.5e0-*w1);
    goto S250;
S110:
    *w = bpser(&a0,&b0,&x0,&eps);
    *w1 = 0.5e0+(0.5e0-*w);
    goto S250;
S120:
    *w1 = bpser(&b0,&a0,&y0,&eps);
    *w = 0.5e0+(0.5e0-*w1);
    goto S250;
S130:
    T2 = 15.0e0*eps;
    *w = bfrac(&a0,&b0,&x0,&y0,&lambda,&T2);
    *w1 = 0.5e0+(0.5e0-*w);
    goto S250;
S140:
    *w1 = bup(&b0,&a0,&y0,&x0,&n,&eps);
    b0 += (double)n;
S150:
    T3 = 15.0e0*eps;
    bgrat(&b0,&a0,&y0,&x0,w1,&T3,&ierr1);
    *w = 0.5e0+(0.5e0-*w1);
    goto S250;
S160:
    n = b0;
    b0 -= (double)n;
    if(b0 != 0.0e0) goto S170;
    n -= 1;
    b0 = 1.0e0;
S170:
    *w = bup(&b0,&a0,&y0,&x0,&n,&eps);
    if(x0 > 0.7e0) goto S180;
    *w += bpser(&a0,&b0,&x0,&eps);
    *w1 = 0.5e0+(0.5e0-*w);
    goto S250;
S180:
    if(a0 > 15.0e0) goto S190;
    n = 20;
    *w += bup(&a0,&b0,&x0,&y0,&n,&eps);
    a0 += (double)n;
S190:
    T4 = 15.0e0*eps;
    bgrat(&a0,&b0,&x0,&y0,w,&T4,&ierr1);
    *w1 = 0.5e0+(0.5e0-*w);
    goto S250;
S200:
    T5 = 100.0e0*eps;
    *w = basym(&a0,&b0,&lambda,&T5);
    *w1 = 0.5e0+(0.5e0-*w);
    goto S250;
S210:
/*
               TERMINATION OF THE PROCEDURE
*/
    if(*a == 0.0e0) goto S320;
S220:
    *w = 0.0e0;
    *w1 = 1.0e0;
    return;
S230:
    if(*b == 0.0e0) goto S330;
S240:
    *w = 1.0e0;
    *w1 = 0.0e0;
    return;
S250:
    if(ind == 0) return;
    t = *w;
    *w = *w1;
    *w1 = t;
    return;
S260:
/*
           PROCEDURE FOR A AND B .LT. 1.E-3*EPS
*/
    *w = *b/(*a+*b);
    *w1 = *a/(*a+*b);
    return;
S270:
/*
                       ERROR RETURN
*/
    *ierr = 1;
    return;
S280:
    *ierr = 2;
    return;
S290:
    *ierr = 3;
    return;
S300:
    *ierr = 4;
    return;
S310:
    *ierr = 5;
    return;
S320:
    *ierr = 6;
    return;
S330:
    *ierr = 7;
    return;
}
double brcmp1(int *mu,double *a,double *b,double *x,double *y)
/*
-----------------------------------------------------------------------
          EVALUATION OF  EXP(MU) * (X**A*Y**B/BETA(A,B))
-----------------------------------------------------------------------
*/
{
static double Const = .398942280401433e0;
static double brcmp1,a0,apb,b0,c,e,h,lambda,lnx,lny,t,u,v,x0,y0,z;
static int i,n;
/*
-----------------
     CONST = 1/SQRT(2*PI)
-----------------
*/
static double T1,T2,T3,T4;
/*
     ..
     .. Executable Statements ..
*/
    a0 = fifdmin1(*a,*b);
    if(a0 >= 8.0e0) goto S130;
    if(*x > 0.375e0) goto S10;
    lnx = log(*x);
    T1 = -*x;
    lny = alnrel(&T1);
    goto S30;
S10:
    if(*y > 0.375e0) goto S20;
    T2 = -*y;
    lnx = alnrel(&T2);
    lny = log(*y);
    goto S30;
S20:
    lnx = log(*x);
    lny = log(*y);
S30:
    z = *a*lnx+*b*lny;
    if(a0 < 1.0e0) goto S40;
    z -= betaln(a,b);
    brcmp1 = esum(mu,&z);
    return brcmp1;
S40:
/*
-----------------------------------------------------------------------
              PROCEDURE FOR A .LT. 1 OR B .LT. 1
-----------------------------------------------------------------------
*/
    b0 = fifdmax1(*a,*b);
    if(b0 >= 8.0e0) goto S120;
    if(b0 > 1.0e0) goto S70;
/*
                   ALGORITHM FOR B0 .LE. 1
*/
    brcmp1 = esum(mu,&z);
    if(brcmp1 == 0.0e0) return brcmp1;
    apb = *a+*b;
    if(apb > 1.0e0) goto S50;
    z = 1.0e0+gam1(&apb);
    goto S60;
S50:
    u = *a+*b-1.e0;
    z = (1.0e0+gam1(&u))/apb;
S60:
    c = (1.0e0+gam1(a))*(1.0e0+gam1(b))/z;
    brcmp1 = brcmp1*(a0*c)/(1.0e0+a0/b0);
    return brcmp1;
S70:
/*
                ALGORITHM FOR 1 .LT. B0 .LT. 8
*/
    u = gamln1(&a0);
    n = b0-1.0e0;
    if(n < 1) goto S90;
    c = 1.0e0;
    for(i=1; i<=n; i++) {
        b0 -= 1.0e0;
        c *= (b0/(a0+b0));
    }
    u = log(c)+u;
S90:
    z -= u;
    b0 -= 1.0e0;
    apb = a0+b0;
    if(apb > 1.0e0) goto S100;
    t = 1.0e0+gam1(&apb);
    goto S110;
S100:
    u = a0+b0-1.e0;
    t = (1.0e0+gam1(&u))/apb;
S110:
    brcmp1 = a0*esum(mu,&z)*(1.0e0+gam1(&b0))/t;
    return brcmp1;
S120:
/*
                   ALGORITHM FOR B0 .GE. 8
*/
    u = gamln1(&a0)+algdiv(&a0,&b0);
    T3 = z-u;
    brcmp1 = a0*esum(mu,&T3);
    return brcmp1;
S130:
/*
-----------------------------------------------------------------------
              PROCEDURE FOR A .GE. 8 AND B .GE. 8
-----------------------------------------------------------------------
*/
    if(*a > *b) goto S140;
    h = *a/ *b;
    x0 = h/(1.0e0+h);
    y0 = 1.0e0/(1.0e0+h);
    lambda = *a-(*a+*b)**x;
    goto S150;
S140:
    h = *b/ *a;
    x0 = 1.0e0/(1.0e0+h);
    y0 = h/(1.0e0+h);
    lambda = (*a+*b)**y-*b;
S150:
    e = -(lambda/ *a);
    if(fabs(e) > 0.6e0) goto S160;
    u = rlog1(&e);
    goto S170;
S160:
    u = e-log(*x/x0);
S170:
    e = lambda/ *b;
    if(fabs(e) > 0.6e0) goto S180;
    v = rlog1(&e);
    goto S190;
S180:
    v = e-log(*y/y0);
S190:
    T4 = -(*a*u+*b*v);
    z = esum(mu,&T4);
    brcmp1 = Const*sqrt(*b*x0)*z*exp(-bcorr(a,b));
    return brcmp1;
}
double brcomp(double *a,double *b,double *x,double *y)
/*
-----------------------------------------------------------------------
               EVALUATION OF X**A*Y**B/BETA(A,B)