nifticdf.c

来自「DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.」· C语言 代码 · 共 2,624 行 · 第 1/5 页

/***=====================================================================***/
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;
} /* END */

/***=====================================================================***/
double brcomp(double *a,double *b,double *x,double *y)
/*
-----------------------------------------------------------------------
               EVALUATION OF X**A*Y**B/BETA(A,B)
-----------------------------------------------------------------------
*/
{
static double Const = .398942280401433e0;
static double brcomp,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;
/*
     ..
     .. Executable Statements ..
*/
    brcomp = 0.0e0;
    if(*x == 0.0e0 || *y == 0.0e0) return brcomp;
    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);
    brcomp = exp(z);
    return brcomp;
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
*/
    brcomp = exp(z);
    if(brcomp == 0.0e0) return brcomp;
    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;
    brcomp = brcomp*(a0*c)/(1.0e0+a0/b0);
    return brcomp;
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:
    brcomp = a0*exp(z)*(1.0e0+gam1(&b0))/t;
    return brcomp;
S120:
/*
                   ALGORITHM FOR B0 .GE. 8
*/
    u = gamln1(&a0)+algdiv(&a0,&b0);
    brcomp = a0*exp(z-u);
    return brcomp;
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:
    z = exp(-(*a*u+*b*v));
    brcomp = Const*sqrt(*b*x0)*z*exp(-bcorr(a,b));
    return brcomp;
} /* END */

/***=====================================================================***/
double bup(double *a,double *b,double *x,double *y,int *n,double *eps)
/*
-----------------------------------------------------------------------
     EVALUATION OF IX(A,B) - IX(A+N,B) WHERE N IS A POSITIVE INTEGER.
     EPS IS THE TOLERANCE USED.
-----------------------------------------------------------------------
*/
{
static int K1 = 1;
static int K2 = 0;
static double bup,ap1,apb,d,l,r,t,w;
static int i,k,kp1,mu,nm1;
/*
     ..
     .. Executable Statements ..
*/
/*
          OBTAIN THE SCALING FACTOR EXP(-MU) AND
             EXP(MU)*(X**A*Y**B/BETA(A,B))/A
*/
    apb = *a+*b;
    ap1 = *a+1.0e0;
    mu = 0;
    d = 1.0e0;
    if(*n == 1 || *a < 1.0e0) goto S10;
    if(apb < 1.1e0*ap1) goto S10;
    mu = fabs(exparg(&K1));
    k = exparg(&K2);
    if(k < mu) mu = k;
    t = mu;
    d = exp(-t);
S10:
    bup = brcmp1(&mu,a,b,x,y)/ *a;
    if(*n == 1 || bup == 0.0e0) return bup;
    nm1 = *n-1;
    w = d;
/*
          LET K BE THE INDEX OF THE MAXIMUM TERM
*/
    k = 0;
    if(*b <= 1.0e0) goto S50;
    if(*y > 1.e-4) goto S20;
    k = nm1;
    goto S30;
S20:
    r = (*b-1.0e0)**x/ *y-*a;
    if(r < 1.0e0) goto S50;
    k = t = nm1;
    if(r < t) k = r;
S30:
/*
          ADD THE INCREASING TERMS OF THE SERIES
*/
    for(i=1; i<=k; i++) {
        l = i-1;
        d = (apb+l)/(ap1+l)**x*d;
        w += d;
    }
    if(k == nm1) goto S70;
S50:
/*
          ADD THE REMAINING TERMS OF THE SERIES
*/
    kp1 = k+1;
    for(i=kp1; i<=nm1; i++) {
        l = i-1;
        d = (apb+l)/(ap1+l)**x*d;
        w += d;
        if(d <= *eps*w) goto S70;
    }
S70:
/*
               TERMINATE THE PROCEDURE
*/
    bup *= w;
    return bup;
} /* END */

/***=====================================================================***/
void cdfbet(int *which,double *p,double *q,double *x,double *y,
            double *a,double *b,int *status,double *bound)
/**********************************************************************

      void cdfbet(int *which,double *p,double *q,double *x,double *y,
            double *a,double *b,int *status,double *bound)

               Cumulative Distribution Function
                         BETa Distribution


                              Function


     Calculates any one parameter of the beta distribution given
     values for the others.


                              Arguments


     WHICH --> Integer indicating which of the next four argument
               values is to be calculated from the others.
               Legal range: 1..4
               iwhich = 1 : Calculate P and Q from X,Y,A and B
               iwhich = 2 : Calculate X and Y from P,Q,A and B
               iwhich = 3 : Calculate A from P,Q,X,Y and B
               iwhich = 4 : Calculate B from P,Q,X,Y and A

     P <--> The integral from 0 to X of the chi-square
            distribution.
            Input range: [0, 1].

     Q <--> 1-P.
            Input range: [0, 1].
            P + Q = 1.0.

     X <--> Upper limit of integration of beta density.
            Input range: [0,1].
            Search range: [0,1]

     Y <--> 1-X.
            Input range: [0,1].
            Search range: [0,1]
            X + Y = 1.0.

     A <--> The first parameter of the beta density.
            Input range: (0, +infinity).
            Search range: [1D-300,1D300]

     B <--> The second parameter of the beta density.
            Input range: (0, +infinity).
            Search range: [1D-300,1D300]

     STATUS <-- 0 if calculation completed correctly
               -I if input parameter number I is out of range
                1 if answer appears to be lower than lowest
                  search bound
                2 if answer appears to be higher than greatest
                  search bound
                3 if P + Q .ne. 1
                4 if X + Y .ne. 1

     BOUND <-- Undefined if STATUS is 0

               Bound exceeded by parameter number I if STATUS
               is negative.

               Lower search bound if STATUS is 1.

               Upper search bound if STATUS is 2.


                              Method


     Cumulative distribution function  (P)  is calculated directly by
     code associated with the following reference.

     DiDinato, A. R. and Morris,  A.   H.  Algorithm 708: Significant
     Digit Computation of the Incomplete  Beta  Function Ratios.  ACM
     Trans. Math.  Softw. 18 (1993), 360-373.

     Computation of other parameters involve a seach for a value that
     produces  the desired  value  of P.   The search relies  on  the
     monotinicity of P with the other parameter.


                              Note


     The beta density is proportional to
               t^(A-1) * (1-t)^(B-1)

**********************************************************************/
{
#define tol (1.0e-8)
#define atol (1.0e-50)
#define zero (1.0e-300)
#define inf 1.0e300
#define one 1.0e0
static int K1 = 1;
static double K2 = 0.0e0;
static double K3 = 1.0e0;
static double K8 = 0.5e0;
static double K9 = 5.0e0;
static double fx,xhi,xlo,cum,ccum,xy,pq;
static unsigned long qhi,qleft,qporq;
static double T4,T5,T6,T7,T10,T11,T12,T13,T14,T15;
/*
     ..
     .. Executable Statements ..
*/
/*
     Check arguments
*/
    if(!(*which < 1 || *which > 4)) goto S30;
    if(!(*which < 1)) goto S10;
    *bound = 1.0e0;
    goto S20;
S10:
    *bound = 4.0e0;
S20:
    *status = -1;
    return;
S30:
    if(*which == 1) goto S70;
/*
     P
*/
    if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
    if(!(*p < 0.0e0)) goto S40;
    *bound = 0.0e0;
    goto S50;
S40:
    *bound = 1.0e0;
S50:
    *status = -2;
    return;
S70:
S60:
    if(*which == 1) goto S110;
/*
     Q
*/
    if(!(*q < 0.0e0 || *q > 1.0e0)) goto S100;
    if(!(*q < 0.0e0)) goto S80;
    *bound = 0.0e0;
    goto S90;
S80:
    *bound = 1.0e0;
S90:
    *status = -3;
    return;
S110:
S100:
    if(*which == 2) goto S150;
/*
     X
*/
    if(!(*x < 0.0e0 || *x > 1.0e0)) goto S140;
    if(!(*x < 0.0e0)) goto S120;
    *bound = 0.0e0;
    goto S130;