dcdflib.c

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

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "cdflib.h"
/*
-----------------------------------------------------------------------
 
     COMPUTATION OF LN(GAMMA(B)/GAMMA(A+B)) WHEN B .GE. 8
 
                         --------
 
     IN THIS ALGORITHM, DEL(X) IS THE FUNCTION DEFINED BY
     LN(GAMMA(X)) = (X - 0.5)*LN(X) - X + 0.5*LN(2*PI) + DEL(X).
 
-----------------------------------------------------------------------
*/
double algdiv(double *a,double *b)
{
static double c0 = .833333333333333e-01;
static double c1 = -.277777777760991e-02;
static double c2 = .793650666825390e-03;
static double c3 = -.595202931351870e-03;
static double c4 = .837308034031215e-03;
static double c5 = -.165322962780713e-02;
static double algdiv,c,d,h,s11,s3,s5,s7,s9,t,u,v,w,x,x2,T1;
/*
     ..
     .. Executable Statements ..
*/
    if(*a <= *b) goto S10;
    h = *b/ *a;
    c = 1.0e0/(1.0e0+h);
    x = h/(1.0e0+h);
    d = *a+(*b-0.5e0);
    goto S20;
S10:
    h = *a/ *b;
    c = h/(1.0e0+h);
    x = 1.0e0/(1.0e0+h);
    d = *b+(*a-0.5e0);
S20:
/*
                SET SN = (1 - X**N)/(1 - X)
*/
    x2 = x*x;
    s3 = 1.0e0+(x+x2);
    s5 = 1.0e0+(x+x2*s3);
    s7 = 1.0e0+(x+x2*s5);
    s9 = 1.0e0+(x+x2*s7);
    s11 = 1.0e0+(x+x2*s9);
/*
                SET W = DEL(B) - DEL(A + B)
*/
    t = pow(1.0e0/ *b,2.0);
    w = ((((c5*s11*t+c4*s9)*t+c3*s7)*t+c2*s5)*t+c1*s3)*t+c0;
    w *= (c/ *b);
/*
                    COMBINE THE RESULTS
*/
    T1 = *a/ *b;
    u = d*alnrel(&T1);
    v = *a*(log(*b)-1.0e0);
    if(u <= v) goto S30;
    algdiv = w-v-u;
    return algdiv;
S30:
    algdiv = w-u-v;
    return algdiv;
}
double alngam(double *x)
/*
**********************************************************************
 
     double alngam(double *x)
                 double precision LN of the GAMma function
 
 
                              Function
 
 
     Returns the natural logarithm of GAMMA(X).
 
 
                              Arguments
 
 
     X --> value at which scaled log gamma is to be returned
                    X is DOUBLE PRECISION
 
 
                              Method
 
 
     If X .le. 6.0, then use recursion to get X below 3
     then apply rational approximation number 5236 of
     Hart et al, Computer Approximations, John Wiley and
     Sons, NY, 1968.
 
     If X .gt. 6.0, then use recursion to get X to at least 12 and
     then use formula 5423 of the same source.
 
**********************************************************************
*/
{
#define hln2pi 0.91893853320467274178e0
static double coef[5] = {
    0.83333333333333023564e-1,-0.27777777768818808e-2,0.79365006754279e-3,
    -0.594997310889e-3,0.8065880899e-3
};
static double scoefd[4] = {
    0.62003838007126989331e2,0.9822521104713994894e1,-0.8906016659497461257e1,
    0.1000000000000000000e1
};
static double scoefn[9] = {
    0.62003838007127258804e2,0.36036772530024836321e2,0.20782472531792126786e2,
    0.6338067999387272343e1,0.215994312846059073e1,0.3980671310203570498e0,
    0.1093115956710439502e0,0.92381945590275995e-2,0.29737866448101651e-2
};
static int K1 = 9;
static int K3 = 4;
static int K5 = 5;
static double alngam,offset,prod,xx;
static int i,n;
static double T2,T4,T6;
/*
     ..
     .. Executable Statements ..
*/
    if(!(*x <= 6.0e0)) goto S70;
    prod = 1.0e0;
    xx = *x;
    if(!(*x > 3.0e0)) goto S30;
S10:
    if(!(xx > 3.0e0)) goto S20;
    xx -= 1.0e0;
    prod *= xx;
    goto S10;
S30:
S20:
    if(!(*x < 2.0e0)) goto S60;
S40:
    if(!(xx < 2.0e0)) goto S50;
    prod /= xx;
    xx += 1.0e0;
    goto S40;
S60:
S50:
    T2 = xx-2.0e0;
    T4 = xx-2.0e0;
    alngam = devlpl(scoefn,&K1,&T2)/devlpl(scoefd,&K3,&T4);
/*
     COMPUTE RATIONAL APPROXIMATION TO GAMMA(X)
*/
    alngam *= prod;
    alngam = log(alngam);
    goto S110;
S70:
    offset = hln2pi;
/*
     IF NECESSARY MAKE X AT LEAST 12 AND CARRY CORRECTION IN OFFSET
*/
    n = fifidint(12.0e0-*x);
    if(!(n > 0)) goto S90;
    prod = 1.0e0;
    for(i=1; i<=n; i++) prod *= (*x+(double)(i-1));
    offset -= log(prod);
    xx = *x+(double)n;
    goto S100;
S90:
    xx = *x;
S100:
/*
     COMPUTE POWER SERIES
*/
    T6 = 1.0e0/pow(xx,2.0);
    alngam = devlpl(coef,&K5,&T6)/xx;
    alngam += (offset+(xx-0.5e0)*log(xx)-xx);
S110:
    return alngam;
#undef hln2pi
}
double alnrel(double *a)
/*
-----------------------------------------------------------------------
            EVALUATION OF THE FUNCTION LN(1 + A)
-----------------------------------------------------------------------
*/
{
static double p1 = -.129418923021993e+01;
static double p2 = .405303492862024e+00;
static double p3 = -.178874546012214e-01;
static double q1 = -.162752256355323e+01;
static double q2 = .747811014037616e+00;
static double q3 = -.845104217945565e-01;
static double alnrel,t,t2,w,x;
/*
     ..
     .. Executable Statements ..
*/
    if(fabs(*a) > 0.375e0) goto S10;
    t = *a/(*a+2.0e0);
    t2 = t*t;
    w = (((p3*t2+p2)*t2+p1)*t2+1.0e0)/(((q3*t2+q2)*t2+q1)*t2+1.0e0);
    alnrel = 2.0e0*t*w;
    return alnrel;
S10:
    x = 1.e0+*a;
    alnrel = log(x);
    return alnrel;
}
double apser(double *a,double *b,double *x,double *eps)
/*
-----------------------------------------------------------------------
     APSER YIELDS THE INCOMPLETE BETA RATIO I(SUB(1-X))(B,A) FOR
     A .LE. MIN(EPS,EPS*B), B*X .LE. 1, AND X .LE. 0.5. USED WHEN
     A IS VERY SMALL. USE ONLY IF ABOVE INEQUALITIES ARE SATISFIED.
-----------------------------------------------------------------------
*/
{
static double g = .577215664901533e0;
static double apser,aj,bx,c,j,s,t,tol;
/*
     ..
     .. Executable Statements ..
*/
    bx = *b**x;
    t = *x-bx;
    if(*b**eps > 2.e-2) goto S10;
    c = log(*x)+psi(b)+g+t;
    goto S20;
S10:
    c = log(bx)+g+t;
S20:
    tol = 5.0e0**eps*fabs(c);
    j = 1.0e0;
    s = 0.0e0;
S30:
    j += 1.0e0;
    t *= (*x-bx/j);
    aj = t/j;
    s += aj;
    if(fabs(aj) > tol) goto S30;
    apser = -(*a*(c+s));
    return apser;
}
double basym(double *a,double *b,double *lambda,double *eps)
/*
-----------------------------------------------------------------------
     ASYMPTOTIC EXPANSION FOR IX(A,B) FOR LARGE A AND B.
     LAMBDA = (A + B)*Y - B  AND EPS IS THE TOLERANCE USED.
     IT IS ASSUMED THAT LAMBDA IS NONNEGATIVE AND THAT
     A AND B ARE GREATER THAN OR EQUAL TO 15.
-----------------------------------------------------------------------
*/
{
static double e0 = 1.12837916709551e0;
static double e1 = .353553390593274e0;
static int num = 20;
/*
------------------------
     ****** NUM IS THE MAXIMUM VALUE THAT N CAN TAKE IN THE DO LOOP
            ENDING AT STATEMENT 50. IT IS REQUIRED THAT NUM BE EVEN.
            THE ARRAYS A0, B0, C, D HAVE DIMENSION NUM + 1.
------------------------
     E0 = 2/SQRT(PI)
     E1 = 2**(-3/2)
------------------------
*/
static int K3 = 1;
static double basym,bsum,dsum,f,h,h2,hn,j0,j1,r,r0,r1,s,sum,t,t0,t1,u,w,w0,z,z0,
    z2,zn,znm1;
static int i,im1,imj,j,m,mm1,mmj,n,np1;
static double a0[21],b0[21],c[21],d[21],T1,T2;
/*
     ..
     .. Executable Statements ..
*/
    basym = 0.0e0;
    if(*a >= *b) goto S10;
    h = *a/ *b;
    r0 = 1.0e0/(1.0e0+h);
    r1 = (*b-*a)/ *b;
    w0 = 1.0e0/sqrt(*a*(1.0e0+h));
    goto S20;
S10:
    h = *b/ *a;
    r0 = 1.0e0/(1.0e0+h);
    r1 = (*b-*a)/ *a;
    w0 = 1.0e0/sqrt(*b*(1.0e0+h));
S20:
    T1 = -(*lambda/ *a);
    T2 = *lambda/ *b;
    f = *a*rlog1(&T1)+*b*rlog1(&T2);
    t = exp(-f);
    if(t == 0.0e0) return basym;
    z0 = sqrt(f);
    z = 0.5e0*(z0/e1);
    z2 = f+f;
    a0[0] = 2.0e0/3.0e0*r1;
    c[0] = -(0.5e0*a0[0]);
    d[0] = -c[0];
    j0 = 0.5e0/e0*erfc1(&K3,&z0);
    j1 = e1;
    sum = j0+d[0]*w0*j1;
    s = 1.0e0;
    h2 = h*h;
    hn = 1.0e0;
    w = w0;
    znm1 = z;
    zn = z2;
    for(n=2; n<=num; n+=2) {
        hn = h2*hn;
        a0[n-1] = 2.0e0*r0*(1.0e0+h*hn)/((double)n+2.0e0);
        np1 = n+1;
        s += hn;
        a0[np1-1] = 2.0e0*r1*s/((double)n+3.0e0);
        for(i=n; i<=np1; i++) {
            r = -(0.5e0*((double)i+1.0e0));
            b0[0] = r*a0[0];
            for(m=2; m<=i; m++) {
                bsum = 0.0e0;
                mm1 = m-1;
                for(j=1; j<=mm1; j++) {
                    mmj = m-j;
                    bsum += (((double)j*r-(double)mmj)*a0[j-1]*b0[mmj-1]);
                }
                b0[m-1] = r*a0[m-1]+bsum/(double)m;
            }
            c[i-1] = b0[i-1]/((double)i+1.0e0);
            dsum = 0.0e0;
            im1 = i-1;
            for(j=1; j<=im1; j++) {
                imj = i-j;
                dsum += (d[imj-1]*c[j-1]);
            }
            d[i-1] = -(dsum+c[i-1]);
        }
        j0 = e1*znm1+((double)n-1.0e0)*j0;
        j1 = e1*zn+(double)n*j1;
        znm1 = z2*znm1;
        zn = z2*zn;
        w = w0*w;
        t0 = d[n-1]*w*j0;
        w = w0*w;
        t1 = d[np1-1]*w*j1;
        sum += (t0+t1);
        if(fabs(t0)+fabs(t1) <= *eps*sum) goto S80;
    }
S80:
    u = exp(-bcorr(a,b));
    basym = e0*t*u*sum;
    return basym;
}
double bcorr(double *a0,double *b0)
/*
-----------------------------------------------------------------------
 
     EVALUATION OF  DEL(A0) + DEL(B0) - DEL(A0 + B0)  WHERE
     LN(GAMMA(A)) = (A - 0.5)*LN(A) - A + 0.5*LN(2*PI) + DEL(A).
     IT IS ASSUMED THAT A0 .GE. 8 AND B0 .GE. 8.
 
-----------------------------------------------------------------------
*/
{
static double c0 = .833333333333333e-01;
static double c1 = -.277777777760991e-02;
static double c2 = .793650666825390e-03;
static double c3 = -.595202931351870e-03;
static double c4 = .837308034031215e-03;
static double c5 = -.165322962780713e-02;
static double bcorr,a,b,c,h,s11,s3,s5,s7,s9,t,w,x,x2;
/*
     ..
     .. Executable Statements ..
*/
    a = fifdmin1(*a0,*b0);
    b = fifdmax1(*a0,*b0);
    h = a/b;
    c = h/(1.0e0+h);
    x = 1.0e0/(1.0e0+h);
    x2 = x*x;
/*
                SET SN = (1 - X**N)/(1 - X)
*/
    s3 = 1.0e0+(x+x2);
    s5 = 1.0e0+(x+x2*s3);
    s7 = 1.0e0+(x+x2*s5);
    s9 = 1.0e0+(x+x2*s7);
    s11 = 1.0e0+(x+x2*s9);
/*
                SET W = DEL(B) - DEL(A + B)
*/
    t = pow(1.0e0/b,2.0);
    w = ((((c5*s11*t+c4*s9)*t+c3*s7)*t+c2*s5)*t+c1*s3)*t+c0;
    w *= (c/b);
/*
                   COMPUTE  DEL(A) + W
*/
    t = pow(1.0e0/a,2.0);
    bcorr = (((((c5*t+c4)*t+c3)*t+c2)*t+c1)*t+c0)/a+w;
    return bcorr;
}
double betaln(double *a0,double *b0)
/*
-----------------------------------------------------------------------
     EVALUATION OF THE LOGARITHM OF THE BETA FUNCTION
-----------------------------------------------------------------------
     E = 0.5*LN(2*PI)
--------------------------
*/
{
static double e = .918938533204673e0;
static double betaln,a,b,c,h,u,v,w,z;
static int i,n;
static double T1;
/*
     ..
     .. Executable Statements ..
*/
    a = fifdmin1(*a0,*b0);
    b = fifdmax1(*a0,*b0);
    if(a >= 8.0e0) goto S100;
    if(a >= 1.0e0) goto S20;
/*
-----------------------------------------------------------------------
                   PROCEDURE WHEN A .LT. 1
-----------------------------------------------------------------------
*/
    if(b >= 8.0e0) goto S10;
    T1 = a+b;
    betaln = gamln(&a)+(gamln(&b)-gamln(&T1));
    return betaln;
S10:
    betaln = gamln(&a)+algdiv(&a,&b);
    return betaln;
S20:
/*
-----------------------------------------------------------------------
                PROCEDURE WHEN 1 .LE. A .LT. 8
-----------------------------------------------------------------------
*/
    if(a > 2.0e0) goto S40;
    if(b > 2.0e0) goto S30;
    betaln = gamln(&a)+gamln(&b)-gsumln(&a,&b);
    return betaln;
S30:
    w = 0.0e0;
    if(b < 8.0e0) goto S60;
    betaln = gamln(&a)+algdiv(&a,&b);
    return betaln;
S40:
/*
                REDUCTION OF A WHEN B .LE. 1000
*/
    if(b > 1000.0e0) goto S80;
    n = a-1.0e0;
    w = 1.0e0;
    for(i=1; i<=n; i++) {
        a -= 1.0e0;
        h = a/b;
        w *= (h/(1.0e0+h));
    }
    w = log(w);
    if(b < 8.0e0) goto S60;
    betaln = w+gamln(&a)+algdiv(&a,&b);
    return betaln;
S60:
/*
                 REDUCTION OF B WHEN B .LT. 8
*/
    n = b-1.0e0;
    z = 1.0e0;
    for(i=1; i<=n; i++) {
        b -= 1.0e0;
        z *= (b/(a+b));
    }
    betaln = w+log(z)+(gamln(&a)+(gamln(&b)-gsumln(&a,&b)));
    return betaln;
S80:
/*
                REDUCTION OF A WHEN B .GT. 1000
*/
    n = a-1.0e0;
    w = 1.0e0;
    for(i=1; i<=n; i++) {
        a -= 1.0e0;
        w *= (a/(1.0e0+a/b));
    }
    betaln = log(w)-(double)n*log(b)+(gamln(&a)+algdiv(&a,&b));
    return betaln;
S100:
/*
-----------------------------------------------------------------------
                   PROCEDURE WHEN A .GE. 8
-----------------------------------------------------------------------
*/
    w = bcorr(&a,&b);
    h = a/b;
    c = h/(1.0e0+h);
    u = -((a-0.5e0)*log(c));
    v = b*alnrel(&h);
    if(u <= v) goto S110;
    betaln = -(0.5e0*log(b))+e+w-v-u;
    return betaln;
S110:
    betaln = -(0.5e0*log(b))+e+w-u-v;
    return betaln;
}
double bfrac(double *a,double *b,double *x,double *y,double *lambda,
	     double *eps)
/*
-----------------------------------------------------------------------
     CONTINUED FRACTION EXPANSION FOR IX(A,B) WHEN A,B .GT. 1.
     IT IS ASSUMED THAT  LAMBDA = (A + B)*Y - B.
-----------------------------------------------------------------------
*/
{
static double bfrac,alpha,an,anp1,beta,bn,bnp1,c,c0,c1,e,n,p,r,r0,s,t,w,yp1;
/*
     ..
     .. Executable Statements ..
*/
    bfrac = brcomp(a,b,x,y);
    if(bfrac == 0.0e0) return bfrac;
    c = 1.0e0+*lambda;
    c0 = *b/ *a;
    c1 = 1.0e0+1.0e0/ *a;
    yp1 = *y+1.0e0;
    n = 0.0e0;
    p = 1.0e0;
    s = *a+1.0e0;
    an = 0.0e0;
    bn = anp1 = 1.0e0;
    bnp1 = c/c1;
    r = c1/c;
S10:
/*
        CONTINUED FRACTION CALCULATION
*/
    n += 1.0e0;
    t = n/ *a;
    w = n*(*b-n)**x;
    e = *a/s;
    alpha = p*(p+c0)*e*e*(w**x);
    e = (1.0e0+t)/(c1+t+t);
    beta = n+w/s+e*(c+n*yp1);
    p = 1.0e0+t;
    s += 2.0e0;
/*
        UPDATE AN, BN, ANP1, AND BNP1
*/
    t = alpha*an+beta*anp1;
    an = anp1;
    anp1 = t;
    t = alpha*bn+beta*bnp1;
    bn = bnp1;
    bnp1 = t;
    r0 = r;
    r = anp1/bnp1;
    if(fabs(r-r0) <= *eps*r) goto S20;
/*
        RESCALE AN, BN, ANP1, AND BNP1
*/
    an /= bnp1;
    bn /= bnp1;
    anp1 = r;
    bnp1 = 1.0e0;