2 : P + 1            - mean vector               P+2 : P*(P+3)/2 + 1  - upper half of cholesky                                       decomposition of cov matrix     x    <-- Vector deviate generated.     work <--> Scratch array                              Method     1) Generate P independent standard normal deviates - Ei ~ N(0,1)     2) Using Cholesky decomposition find A s.t. trans(A)*A = COVM     3) trans(A)E + MEANV ~ N(MEANV,COVM)***********************************************************************/{static long i,icount,j,p,D1,D2,D3,D4;static float ae;    p = (long) (*parm);/*     Generate P independent normal deviates - WORK ~ N(0,1)*/    for(i=1; i<=p; i++) *(work+i-1) = snorm();    for(i=1,D3=1,D4=(p-i+D3)/D3; D4>0; D4--,i+=D3) {/*     PARM (P+2 : P*(P+3)/2 + 1) contains A, the Cholesky      decomposition of the desired covariance matrix.          trans(A)(1,1) = PARM(P+2)          trans(A)(2,1) = PARM(P+3)          trans(A)(2,2) = PARM(P+2+P)          trans(A)(3,1) = PARM(P+4)          trans(A)(3,2) = PARM(P+3+P)          trans(A)(3,3) = PARM(P+2-1+2P)  ...     trans(A)*WORK + MEANV ~ N(MEANV,COVM)*/        icount = 0;        ae = 0.0;        for(j=1,D1=1,D2=(i-j+D1)/D1; D2>0; D2--,j+=D1) {            icount += (j-1);            ae += (*(parm+i+(j-1)*p-icount+p)**(work+j-1));        }        *(x+i-1) = ae+*(parm+i);    }}void genmul(long n,float *p,long ncat,long *ix)/***********************************************************************      void genmul(int n,float *p,int ncat,int *ix)     GENerate an observation from the MULtinomial distribution                              Arguments     N --> Number of events that will be classified into one of           the categories 1..NCAT     P --> Vector of probabilities.  P(i) is the probability that           an event will be classified into category i.  Thus, P(i)           must be [0,1]. Only the first NCAT-1 P(i) must be defined           since P(NCAT) is 1.0 minus the sum of the first           NCAT-1 P(i).     NCAT --> Number of categories.  Length of P and IX.     IX <-- Observation from multinomial distribution.  All IX(i)            will be nonnegative and their sum will be N.                              Method     Algorithm from page 559 of      Devroye, Luc      Non-Uniform Random Variate Generation.  Springer-Verlag,     New York, 1986. ***********************************************************************/{static float prob,ptot,sum;static long i,icat,ntot;    if(n < 0) ftnstop("N < 0 in GENMUL");    if(ncat <= 1) ftnstop("NCAT <= 1 in GENMUL");    ptot = 0.0F;    for(i=0; i<ncat-1; i++) {        if(*(p+i) < 0.0F) ftnstop("Some P(i) < 0 in GENMUL");        if(*(p+i) > 1.0F) ftnstop("Some P(i) > 1 in GENMUL");        ptot += *(p+i);    }    if(ptot > 0.99999F) ftnstop("Sum of P(i) > 1 in GENMUL");/*     Initialize variables*/    ntot = n;    sum = 1.0F;    for(i=0; i<ncat; i++) ix[i] = 0;/*     Generate the observation*/    for(icat=0; icat<ncat-1; icat++) {        prob = *(p+icat)/sum;        *(ix+icat) = ignbin(ntot,prob);        ntot -= *(ix+icat);	if(ntot <= 0) return;        sum -= *(p+icat);    }    *(ix+ncat-1) = ntot;/*     Finished*/    return;}float gennch(float df,float xnonc)/***********************************************************************     float gennch(float df,float xnonc)           Generate random value of Noncentral CHIsquare variable                              Function     Generates random deviate  from the  distribution  of a  noncentral     chisquare with DF degrees  of freedom and noncentrality  parameter     xnonc.                              Arguments     df --> Degrees of freedom of the chisquare            (Must be >= 1.0)     xnonc --> Noncentrality parameter of the chisquare               (Must be >= 0.0)                              Method     Uses fact that  noncentral chisquare  is  the  sum of a  chisquare     deviate with DF-1  degrees of freedom plus the  square of a normal     deviate with mean XNONC and standard deviation 1.***********************************************************************/{static float gennch;    if(!(df < 1.0 || xnonc < 0.0)) goto S10;    fputs("DF < 1 or XNONC < 0 in GENNCH - ABORT\n",stderr);    fprintf(stderr,"Value of DF: %16.6E Value of XNONC: %16.6E\n",df,xnonc);    exit(1);/* JJV changed code to call SGAMMA, SNORM directly */S10:    if(df >= 1.000001) goto S20;/* * JJV case df == 1.0 * gennch = pow(gennor(sqrt(xnonc),1.0),2.0); <- OLD */    gennch = pow(snorm()+sqrt(xnonc),2.0);    goto S30;S20:/* * JJV case df > 1.0 * gennch = genchi(df-1.0)+pow(gennor(sqrt(xnonc),1.0),2.0); <- OLD */    gennch = 2.0*sgamma((df-1.0)/2.0)+pow(snorm()+sqrt(xnonc),2.0);S30:    return gennch;}float gennf(float dfn,float dfd,float xnonc)/***********************************************************************     float gennf(float dfn,float dfd,float xnonc)           GENerate random deviate from the Noncentral F distribution                              Function     Generates a random deviate from the  noncentral F (variance ratio)     distribution with DFN degrees of freedom in the numerator, and DFD     degrees of freedom in the denominator, and noncentrality parameter     XNONC.                              Arguments     dfn --> Numerator degrees of freedom             (Must be >= 1.0)     dfd --> Denominator degrees of freedom             (Must be positive)     xnonc --> Noncentrality parameter               (Must be nonnegative)                              Method     Directly generates ratio of noncentral numerator chisquare variate     to central denominator chisquare variate.***********************************************************************/{static float gennf,xden,xnum;static long qcond;    /* JJV changed qcond, error message to allow dfn == 1.0 */    qcond = dfn < 1.0 || dfd <= 0.0 || xnonc < 0.0;    if(!qcond) goto S10;    fputs("In GENNF - Either (1) Numerator DF < 1.0 or\n",stderr);    fputs(" (2) Denominator DF <= 0.0 or\n",stderr);    fputs(" (3) Noncentrality parameter < 0.0\n",stderr);    fprintf(stderr,      "DFN value: %16.6E DFD value: %16.6E XNONC value: \n%16.6E\n",dfn,dfd,      xnonc);    exit(1);S10:/* * JJV changed the code to call SGAMMA and SNORM directly * GENNF = ( GENNCH( DFN, XNONC ) / DFN ) / ( GENCHI( DFD ) / DFD ) * xnum = gennch(dfn,xnonc)/dfn; <- OLD * xden = genchi(dfd)/dfd; <- OLD */    if(dfn >= 1.000001) goto S20;/* JJV case dfn == 1.0, dfn is counted as exactly 1.0 */    xnum = pow(snorm()+sqrt(xnonc),2.0);    goto S30;S20:/* JJV case df > 1.0 */    xnum = (2.0*sgamma((dfn-1.0)/2.0)+pow(snorm()+sqrt(xnonc),2.0))/dfn;S30:    xden = 2.0*sgamma(dfd/2.0)/dfd;/* * JJV changed constant to prevent underflow at compile time. *   if(!(xden <= 9.999999999998E-39*xnum)) goto S40; */    if(!(xden <= 1.0E-37*xnum)) goto S40;    fputs(" GENNF - generated numbers would cause overflow\n",stderr);    fprintf(stderr," Numerator %16.6E Denominator %16.6E\n",xnum,xden);/* * JJV changed next 2 lines to reflect constant change above in the * JJV truncated value returned. *   fputs(" GENNF returning 1.0E38\n",stderr); *   gennf = 1.0E38; */    fputs(" GENNF returning 1.0E37\n",stderr);    gennf = 1.0E37;    goto S50;S40:    gennf = xnum/xden;S50:    return gennf;}float gennor(float av,float sd)/***********************************************************************     float gennor(float av,float sd)         GENerate random deviate from a NORmal distribution                              Function     Generates a single random deviate from a normal distribution     with mean, AV, and standard deviation, SD.                              Arguments     av --> Mean of the normal distribution.     sd --> Standard deviation of the normal distribution.     JJV    (sd >= 0)                              Method     Renames SNORM from TOMS as slightly modified by BWB to use RANF     instead of SUNIF.     For details see:               Ahrens, J.H. and Dieter, U.               Extensions of Forsythe's Method for Random               Sampling from the Normal Distribution.               Math. Comput., 27,124 (Oct. 1973), 927 - 937.***********************************************************************/{static float gennor;/* JJV added argument checker */    if(sd >= 0.0) goto S10;    fputs(" SD < 0 in GENNOR - ABORT\n",stderr);    fprintf(stderr," Value of SD: %16.6E\n",sd);    exit(1);S10:    gennor = sd*snorm()+av;    return gennor;}void genprm(long *iarray,int larray)/***********************************************************************    void genprm(long *iarray,int larray)               GENerate random PeRMutation of iarray                              Arguments     iarray <--> On output IARRAY is a random permutation of its                 value on input     larray <--> Length of IARRAY***********************************************************************/{static long i,itmp,iwhich,D1,D2;    for(i=1,D1=1,D2=(larray-i+D1)/D1; D2>0; D2--,i+=D1) {        iwhich = ignuin(i,larray);        itmp = *(iarray+iwhich-1);        *(iarray+iwhich-1) = *(iarray+i-1);        *(iarray+i-1) = itmp;    }}float genunf(float low,float high)/***********************************************************************     float genunf(float low,float high)               GeNerate Uniform Real between LOW and HIGH                              Function     Generates a real uniformly distributed between LOW and HIGH.                              Arguments     low --> Low bound (exclusive) on real value to be generated     high --> High bound (exclusive) on real value to be generated***********************************************************************/{static float genunf;    if(!(low > high)) goto S10;    fprintf(stderr,"LOW > HIGH in GENUNF: LOW %16.6E HIGH: %16.6E\n",low,high);    fputs("Abort\n",stderr);    exit(1);S10:    genunf = low+(high-low)*ranf();    return genunf;}void gscgn(long getset,long *g)/***********************************************************************     void gscgn(long getset,long *g)                         Get/Set GeNerator     Gets or returns in G the number of the current generator                              Arguments     getset --> 0 Get                1 Set     g <-- Number of the current random number generator (1..32)***********************************************************************/{#define numg 32Lstatic long curntg = 1;    if(getset == 0) *g = curntg;    else  {        if(*g < 0 || *g > numg) {            fputs(" Generator number out of range in GSCGN\n",stderr);            exit(0);        }        curntg = *g;    }#undef numg}void gsrgs(long getset,long *qvalue)/***********************************************************************     void gsrgs(long getset,long *qvalue)               Get/Set Random Generators Set     Gets or sets whether random generators set (initialized).     Initially (data statement) state is not set     If getset is 1 state is set to qvalue     If getset is 0 state returned in qvalue***********************************************************************/{static long qinit = 0;    if(getset == 0) *qvalue = qinit;    else qinit = *qvalue;}void gssst(long getset,long *qset)/***********************************************************************     void gssst(long getset,long *qset)          Get or Set whether Seed is Set     Initialize to Seed not Set     If getset is 1 sets state to Seed Set     If getset is 0 returns T in qset if Seed Set     Else returns F in qset***********************************************************************/{static long qstate = 0;    if(getset != 0) qstate = 1;    else  *qset = qstate;}long ignbin(long n,float pp)/***********************************************************************     long ignbin(long n,float pp)                    GENerate BINomial random deviate                              Function     Generates a single random deviate from a binomial     distribution whose number of trials is N and whose     probability of an event in each trial is P.                              Arguments     n  --> The number of trials in the binomial distribution            from which a random deviate is to be generated.	    JJV (N >= 0)     pp --> The probability of an event in each trial of the            binomial distribution from which a random deviate            is to be generated.	    JJV (0.0 <= PP <= 1.0)     ignbin <-- A random deviate yielding the number of events                from N independent trials, each of which has                a probability of event P.                              Method     This is algorithm BTPE from:         Kachitvichyanukul, V. and Schmeiser, B. W.         Binomial Random Variate Generation.         Communications of the ACM, 31, 2         (February, 1988) 216.**********************************************************************     SUBROUTINE BTPEC(N,PP,ISEED,JX)     BINOMIAL RANDOM VARIATE GENERATOR     MEAN .LT. 30 -- INVERSE CDF       MEAN .GE. 30 -- ALGORITHM BTPE:  ACCEPTANCE-REJECTION VIA       FOUR REGION COMPOSITION.  THE FOUR REGIONS ARE A TRIANGLE       (SYMMETRIC IN THE CENTER), A PAIR OF PARALLELOGRAMS (ABOVE       THE TRIANGLE), AND EXPONENTIAL LEFT AND RIGHT TAILS.     BTPE REFERS TO BINOMIAL-TRIANGLE-PARALLELOGRAM-EXPONENTIAL.     BTPEC REFERS TO BTPE AND "COMBINED."  THUS BTPE IS THE       RESEARCH AND BTPEC IS THE IMPLEMENTATION OF A COMPLETE       USABLE ALGORITHM.     REFERENCE:  VORATAS KACHITVICHYANUKUL AND BRUCE SCHMEISER,       "BINOMIAL RANDOM VARIATE GENERATION,"       COMMUNICATIONS OF THE ACM, FORTHCOMING     WRITTEN:  SEPTEMBER 1980.       LAST REVISED:  MAY 1985, JULY 1987     REQUIRED SUBPROGRAM:  RAND() -- A UNIFORM (0,1) RANDOM NUMBER                           GENERATOR     ARGUMENTS       N : NUMBER OF BERNOULLI TRIALS            (INPUT)       PP : PROBABILITY OF SUCCESS IN EACH TRIAL (INPUT)       ISEED:  RANDOM NUMBER SEED                (INPUT AND OUTPUT)       JX:  RANDOMLY GENERATED OBSERVATION       (OUTPUT)     VARIABLES       PSAVE: VALUE OF PP FROM THE LAST CALL TO BTPEC       NSAVE: VALUE OF N FROM THE LAST CALL TO BTPEC       XNP:  VALUE OF THE MEAN FROM THE LAST CALL TO BTPEC       P: PROBABILITY USED IN THE GENERATION PHASE OF BTPEC       FFM: TEMPORARY VARIABLE EQUAL TO XNP + P       M:  INTEGER VALUE OF THE CURRENT MODE       FM:  FLOATING POINT VALUE OF THE CURRENT MODE       XNPQ: TEMPORARY VARIABLE USED IN SETUP AND SQUEEZING STEPS       P1:  AREA OF THE TRIANGLE       C:  HEIGHT OF THE PARALLELOGRAMS       XM:  CENTER OF THE TRIANGLE       XL:  LEFT END OF THE TRIANGLE