1:
2: #include "ranlib.h"
3: #include <stdio.h>
4: #include <math.h>
5: #include <stdlib.h>
6: #define abs(x) ((x) >= 0 ? (x) : -(x))
7: #define min(a,b) ((a) <= (b) ? (a) : (b))
8: #define max(a,b) ((a) >= (b) ? (a) : (b))
9:
10: /*
11: **********************************************************************
12: float genbet(float aa,float bb)
13: GeNerate BETa random deviate
14: Function
15: Returns a single random deviate from the beta distribution with
16: parameters A and B. The density of the beta is
17: x^(a-1) * (1-x)^(b-1) / B(a,b) for 0 < x < 1
18: Arguments
19: aa --> First parameter of the beta distribution
20:
21: bb --> Second parameter of the beta distribution
22:
23: Method
24: R. C. H. Cheng
25: Generating Beta Variatew with Nonintegral Shape Parameters
26: Communications of the ACM, 21:317-322 (1978)
27: (Algorithms BB and BC)
28: **********************************************************************
29: */
30: float genbet(float aa,float bb)
31: {
32: #define expmax 89.0
33: #define infnty 1.0E38
34: static float olda = -1.0;
35: static float oldb = -1.0;
36: static float genbet,a,alpha,b,beta,delta,gamma,k1,k2,r,s,t,u1,u2,v,w,y,z;
37: static long qsame;
38:
39: qsame = olda == aa && oldb == bb;
40: if(qsame) goto S20;
41: if(!(aa <= 0.0 || bb <= 0.0)) goto S10;
42: puts(" AA or BB <= 0 in GENBET - Abort!");
43: printf(" AA: %16.6E BB %16.6E\n",aa,bb);
44: exit(1);
45: S10:
46: olda = aa;
47: oldb = bb;
48: S20:
49: if(!(min(aa,bb) > 1.0)) goto S100;
50: /*
51: Alborithm BB
52: Initialize
53: */
54: if(qsame) goto S30;
55: a = min(aa,bb);
56: b = max(aa,bb);
57: alpha = a+b;
58: beta = sqrt((alpha-2.0)/(2.0*a*b-alpha));
59: gamma = a+1.0/beta;
60: S30:
61: S40:
62: u1 = ranf();
63: /*
64: Step 1
65: */
66: u2 = ranf();
67: v = beta*log(u1/(1.0-u1));
68: if(!(v > expmax)) goto S50;
69: w = infnty;
70: goto S60;
71: S50:
72: w = a*exp(v);
73: S60:
74: z = pow(u1,2.0)*u2;
75: r = gamma*v-1.3862944;
76: s = a+r-w;
77: /*
78: Step 2
79: */
80: if(s+2.609438 >= 5.0*z) goto S70;
81: /*
82: Step 3
83: */
84: t = log(z);
85: if(s > t) goto S70;
86: /*
87: Step 4
88: */
89: if(r+alpha*log(alpha/(b+w)) < t) goto S40;
90: S70:
91: /*
92: Step 5
93: */
94: if(!(aa == a)) goto S80;
95: genbet = w/(b+w);
96: goto S90;
97: S80:
98: genbet = b/(b+w);
99: S90:
100: goto S230;
101: S100:
102: /*
103: Algorithm BC
104: Initialize
105: */
106: if(qsame) goto S110;
107: a = max(aa,bb);
108: b = min(aa,bb);
109: alpha = a+b;
110: beta = 1.0/b;
111: delta = 1.0+a-b;
112: k1 = delta*(1.38889E-2+4.16667E-2*b)/(a*beta-0.777778);
113: k2 = 0.25+(0.5+0.25/delta)*b;
114: S110:
115: S120:
116: u1 = ranf();
117: /*
118: Step 1
119: */
120: u2 = ranf();
121: if(u1 >= 0.5) goto S130;
122: /*
123: Step 2
124: */
125: y = u1*u2;
126: z = u1*y;
127: if(0.25*u2+z-y >= k1) goto S120;
128: goto S170;
129: S130:
130: /*
131: Step 3
132: */
133: z = pow(u1,2.0)*u2;
134: if(!(z <= 0.25)) goto S160;
135: v = beta*log(u1/(1.0-u1));
136: if(!(v > expmax)) goto S140;
137: w = infnty;
138: goto S150;
139: S140:
140: w = a*exp(v);
141: S150:
142: goto S200;
143: S160:
144: if(z >= k2) goto S120;
145: S170:
146: /*
147: Step 4
148: Step 5
149: */
150: v = beta*log(u1/(1.0-u1));
151: if(!(v > expmax)) goto S180;
152: w = infnty;
153: goto S190;
154: S180:
155: w = a*exp(v);
156: S190:
157: if(alpha*(log(alpha/(b+w))+v)-1.3862944 < log(z)) goto S120;
158: S200:
159: /*
160: Step 6
161: */
162: if(!(a == aa)) goto S210;
163: genbet = w/(b+w);
164: goto S220;
165: S210:
166: genbet = b/(b+w);
167: S230:
168: S220:
169: return genbet;
170: #undef expmax
171: #undef infnty
172: }
173: float genchi(float df)
174: /*
175: **********************************************************************
176: float genchi(float df)
177: Generate random value of CHIsquare variable
178: Function
179: Generates random deviate from the distribution of a chisquare
180: with DF degrees of freedom random variable.
181: Arguments
182: df --> Degrees of freedom of the chisquare
183: (Must be positive)
184:
185: Method
186: Uses relation between chisquare and gamma.
187: **********************************************************************
188: */
189: {
190: static float genchi;
191:
192: if(!(df <= 0.0)) goto S10;
193: puts("DF <= 0 in GENCHI - ABORT");
194: printf("Value of DF: %16.6E\n",df);
195: exit(1);
196: S10:
197: genchi = 2.0*gengam(1.0,df/2.0);
198: return genchi;
199: }
200: float genexp(float av)
201: /*
202: **********************************************************************
203: float genexp(float av)
204: GENerate EXPonential random deviate
205: Function
206: Generates a single random deviate from an exponential
207: distribution with mean AV.
208: Arguments
209: av --> The mean of the exponential distribution from which
210: a random deviate is to be generated.
211: Method
212: Renames SEXPO from TOMS as slightly modified by BWB to use RANF
213: instead of SUNIF.
214: For details see:
215: Ahrens, J.H. and Dieter, U.
216: Computer Methods for Sampling From the
217: Exponential and Normal Distributions.
218: Comm. ACM, 15,10 (Oct. 1972), 873 - 882.
219: **********************************************************************
220: */
221: {
222: static float genexp;
223:
224: genexp = sexpo()*av;
225: return genexp;
226: }
227: float genf(float dfn,float dfd)
228: /*
229: **********************************************************************
230: float genf(float dfn,float dfd)
231: GENerate random deviate from the F distribution
232: Function
233: Generates a random deviate from the F (variance ratio)
234: distribution with DFN degrees of freedom in the numerator
235: and DFD degrees of freedom in the denominator.
236: Arguments
237: dfn --> Numerator degrees of freedom
238: (Must be positive)
239: dfd --> Denominator degrees of freedom
240: (Must be positive)
241: Method
242: Directly generates ratio of chisquare variates
243: **********************************************************************
244: */
245: {
246: static float genf,xden,xnum;
247:
248: if(!(dfn <= 0.0 || dfd <= 0.0)) goto S10;
249: puts("Degrees of freedom nonpositive in GENF - abort!");
250: printf("DFN value: %16.6EDFD value: %16.6E\n",dfn,dfd);
251: exit(1);
252: S10:
253: xnum = genchi(dfn)/dfn;
254: /*
255: GENF = ( GENCHI( DFN ) / DFN ) / ( GENCHI( DFD ) / DFD )
256: */
257: xden = genchi(dfd)/dfd;
258: if(!(xden <= 9.999999999998E-39*xnum)) goto S20;
259: puts(" GENF - generated numbers would cause overflow");
260: printf(" Numerator %16.6E Denominator %16.6E\n",xnum,xden);
261: puts(" GENF returning 1.0E38");
262: genf = 1.0E38;
263: goto S30;
264: S20:
265: genf = xnum/xden;
266: S30:
267: return genf;
268: }
269: float gengam(float a,float r)
270: /*
271: **********************************************************************
272: float gengam(float a,float r)
273: GENerates random deviates from GAMma distribution
274: Function
275: Generates random deviates from the gamma distribution whose
276: density is
277: (A**R)/Gamma(R) * X**(R-1) * Exp(-A*X)
278: Arguments
279: a --> Location parameter of Gamma distribution
280: r --> Shape parameter of Gamma distribution
281: Method
282: Renames SGAMMA from TOMS as slightly modified by BWB to use RANF
283: instead of SUNIF.
284: For details see:
285: (Case R >= 1.0)
286: Ahrens, J.H. and Dieter, U.
287: Generating Gamma Variates by a
288: Modified Rejection Technique.
289: Comm. ACM, 25,1 (Jan. 1982), 47 - 54.
290: Algorithm GD
291: (Case 0.0 <= R <= 1.0)
292: Ahrens, J.H. and Dieter, U.
293: Computer Methods for Sampling from Gamma,
294: Beta, Poisson and Binomial Distributions.
295: Computing, 12 (1974), 223-246/
296: Adapted algorithm GS.
297: **********************************************************************
298: */
299: {
300: static float gengam;
301:
302: gengam = sgamma(r);
303: gengam /= a;
304: return gengam;
305: }
306: void genmn(float *parm,float *x,float *work)
307: /*
308: **********************************************************************
309: void genmn(float *parm,float *x,float *work)
310: GENerate Multivariate Normal random deviate
311: Arguments
312: parm --> Parameters needed to generate multivariate normal
313: deviates (MEANV and Cholesky decomposition of
314: COVM). Set by a previous call to SETGMN.
315: 1 : 1 - size of deviate, P
316: 2 : P + 1 - mean vector
317: P+2 : P*(P+3)/2 + 1 - upper half of cholesky
318: decomposition of cov matrix
319: x <-- Vector deviate generated.
320: work <--> Scratch array
321: Method
322: 1) Generate P independent standard normal deviates - Ei ~ N(0,1)
323: 2) Using Cholesky decomposition find A s.t. trans(A)*A = COVM
324: 3) trans(A)E + MEANV ~ N(MEANV,COVM)
325: **********************************************************************
326: */
327: {
328: static long i,icount,j,p,D1,D2,D3,D4;
329: static float ae;
330:
331: p = (long) (*parm);
332: /*
333: Generate P independent normal deviates - WORK ~ N(0,1)
334: */
335: for(i=1; i<=p; i++) *(work+i-1) = snorm();
336: for(i=1,D3=1,D4=(p-i+D3)/D3; D4>0; D4--,i+=D3) {
337: /*
338: PARM (P+2 : P*(P+3)/2 + 1) contains A, the Cholesky
339: decomposition of the desired covariance matrix.
340: trans(A)(1,1) = PARM(P+2)
341: trans(A)(2,1) = PARM(P+3)
342: trans(A)(2,2) = PARM(P+2+P)
343: trans(A)(3,1) = PARM(P+4)
344: trans(A)(3,2) = PARM(P+3+P)
345: trans(A)(3,3) = PARM(P+2-1+2P) ...
346: trans(A)*WORK + MEANV ~ N(MEANV,COVM)
347: */
348: icount = 0;
349: ae = 0.0;
350: for(j=1,D1=1,D2=(i-j+D1)/D1; D2>0; D2--,j+=D1) {
351: icount += (j-1);
352: ae += (*(parm+i+(j-1)*p-icount+p)**(work+j-1));
353: }
354: *(x+i-1) = ae+*(parm+i);
355: }
356: }
357: float gennch(float df,float xnonc)
358: /*
359: **********************************************************************
360: float gennch(float df,float xnonc)
361: Generate random value of Noncentral CHIsquare variable
362: Function
363: Generates random deviate from the distribution of a noncentral
364: chisquare with DF degrees of freedom and noncentrality parameter
365: xnonc.
366: Arguments
367: df --> Degrees of freedom of the chisquare
368: (Must be > 1.0)
369: xnonc --> Noncentrality parameter of the chisquare
370: (Must be >= 0.0)
371: Method
372: Uses fact that noncentral chisquare is the sum of a chisquare
373: deviate with DF-1 degrees of freedom plus the square of a normal
374: deviate with mean XNONC and standard deviation 1.
375: **********************************************************************
376: */
377: {
378: static float gennch;
379:
380: if(!(df <= 1.0 || xnonc < 0.0)) goto S10;
381: puts("DF <= 1 or XNONC < 0 in GENNCH - ABORT");
382: printf("Value of DF: %16.6E Value of XNONC%16.6E\n",df,xnonc);
383: exit(1);
384: S10:
385: gennch = genchi(df-1.0)+pow(gennor(sqrt(xnonc),1.0),2.0);
386: return gennch;
387: }
388: float gennf(float dfn,float dfd,float xnonc)
389: /*
390: **********************************************************************
391: float gennf(float dfn,float dfd,float xnonc)
392: GENerate random deviate from the Noncentral F distribution
393: Function
394: Generates a random deviate from the noncentral F (variance ratio)
395: distribution with DFN degrees of freedom in the numerator, and DFD
396: degrees of freedom in the denominator, and noncentrality parameter
397: XNONC.
398: Arguments
399: dfn --> Numerator degrees of freedom
400: (Must be >= 1.0)
401: dfd --> Denominator degrees of freedom
402: (Must be positive)
403: xnonc --> Noncentrality parameter
404: (Must be nonnegative)
405: Method
406: Directly generates ratio of noncentral numerator chisquare variate
407: to central denominator chisquare variate.
408: **********************************************************************
409: */
410: {
411: static float gennf,xden,xnum;
412: static long qcond;
413:
414: qcond = dfn <= 1.0 || dfd <= 0.0 || xnonc < 0.0;
415: if(!qcond) goto S10;
416: puts("In GENNF - Either (1) Numerator DF <= 1.0 or");
417: puts("(2) Denominator DF < 0.0 or ");
418: puts("(3) Noncentrality parameter < 0.0");
419: printf("DFN value: %16.6EDFD value: %16.6EXNONC value: \n%16.6E\n",dfn,dfd,
420: xnonc);
421: exit(1);
422: S10:
423: xnum = gennch(dfn,xnonc)/dfn;
424: /*
425: GENNF = ( GENNCH( DFN, XNONC ) / DFN ) / ( GENCHI( DFD ) / DFD )
426: */
427: xden = genchi(dfd)/dfd;
428: if(!(xden <= 9.999999999998E-39*xnum)) goto S20;
429: puts(" GENNF - generated numbers would cause overflow");
430: printf(" Numerator %16.6E Denominator %16.6E\n",xnum,xden);
431: puts(" GENNF returning 1.0E38");
432: gennf = 1.0E38;
433: goto S30;
434: S20:
435: gennf = xnum/xden;
436: S30:
437: return gennf;
438: }
439: float gennor(float av,float sd)
440: /*
441: **********************************************************************
442: float gennor(float av,float sd)
443: GENerate random deviate from a NORmal distribution
444: Function
445: Generates a single random deviate from a normal distribution
446: with mean, AV, and standard deviation, SD.
447: Arguments
448: av --> Mean of the normal distribution.
449: sd --> Standard deviation of the normal distribution.
450: Method
451: Renames SNORM from TOMS as slightly modified by BWB to use RANF
452: instead of SUNIF.
453: For details see:
454: Ahrens, J.H. and Dieter, U.
455: Extensions of Forsythe's Method for Random
456: Sampling from the Normal Distribution.
457: Math. Comput., 27,124 (Oct. 1973), 927 - 937.
458: **********************************************************************
459: */
460: {
461: static float gennor;
462:
463: gennor = sd*snorm()+av;
464: return gennor;
465: }
466: void genprm(long *iarray,int larray)
467: /*
468: **********************************************************************
469: void genprm(long *iarray,int larray)
470: GENerate random PeRMutation of iarray
471: Arguments
472: iarray <--> On output IARRAY is a random permutation of its
473: value on input
474: larray <--> Length of IARRAY
475: **********************************************************************
476: */
477: {
478: static long i,itmp,iwhich,D1,D2;
479:
480: for(i=1,D1=1,D2=(larray-i+D1)/D1; D2>0; D2--,i+=D1) {
481: iwhich = ignuin(i,larray);
482: itmp = *(iarray+iwhich-1);
483: *(iarray+iwhich-1) = *(iarray+i-1);
484: *(iarray+i-1) = itmp;
485: }
486: }
487: float genunf(float low,float high)
488: /*
489: **********************************************************************
490: float genunf(float low,float high)
491: GeNerate Uniform Real between LOW and HIGH
492: Function
493: Generates a real uniformly distributed between LOW and HIGH.
494: Arguments
495: low --> Low bound (exclusive) on real value to be generated
496: high --> High bound (exclusive) on real value to be generated
497: **********************************************************************
498: */
499: {
500: static float genunf;
501:
502: if(!(low > high)) goto S10;
503: printf("LOW > HIGH in GENUNF: LOW %16.6E HIGH: %16.6E\n",low,high);
504: puts("Abort");
505: exit(1);
506: S10:
507: genunf = low+(high-low)*ranf();
508: return genunf;
509: }
510: void gscgn(long getset,long *g)
511: /*
512: **********************************************************************
513: void gscgn(long getset,long *g)
514: Get/Set GeNerator
515: Gets or returns in G the number of the current generator
516: Arguments
517: getset --> 0 Get
518: 1 Set
519: g <-- Number of the current random number generator (1..32)
520: **********************************************************************
521: */
522: {
523: #define numg 32L
524: static long curntg = 1;
525: if(getset == 0) *g = curntg;
526: else {
527: if(*g < 0 || *g > numg) {
528: puts(" Generator number out of range in GSCGN");
529: exit(0);
530: }
531: curntg = *g;
532: }
533: #undef numg
534: }
535: void gsrgs(long getset,long *qvalue)
536: /*
537: **********************************************************************
538: void gsrgs(long getset,long *qvalue)
539: Get/Set Random Generators Set
540: Gets or sets whether random generators set (initialized).
541: Initially (data statement) state is not set
542: If getset is 1 state is set to qvalue
543: If getset is 0 state returned in qvalue
544: **********************************************************************
545: */
546: {
547: static long qinit = 0;
548:
549: if(getset == 0) *qvalue = qinit;
550: else qinit = *qvalue;
551: }
552: void gssst(long getset,long *qset)
553: /*
554: **********************************************************************
555: void gssst(long getset,long *qset)
556: Get or Set whether Seed is Set
557: Initialize to Seed not Set
558: If getset is 1 sets state to Seed Set
559: If getset is 0 returns T in qset if Seed Set
560: Else returns F in qset
561: **********************************************************************
562: */
563: {
564: static long qstate = 0;
565: if(getset != 0) qstate = 1;
566: else *qset = qstate;
567: }
568: long ignbin(long n,float pp)
569: /*
570: **********************************************************************
571: long ignbin(long n,float pp)
572: GENerate BINomial random deviate
573: Function
574: Generates a single random deviate from a binomial
575: distribution whose number of trials is N and whose
576: probability of an event in each trial is P.
577: Arguments
578: n --> The number of trials in the binomial distribution
579: from which a random deviate is to be generated.
580: p --> The probability of an event in each trial of the
581: binomial distribution from which a random deviate
582: is to be generated.
583: ignbin <-- A random deviate yielding the number of events
584: from N independent trials, each of which has
585: a probability of event P.
586: Method
587: This is algorithm BTPE from:
588: Kachitvichyanukul, V. and Schmeiser, B. W.
589: Binomial Random Variate Generation.
590: Communications of the ACM, 31, 2
591: (February, 1988) 216.
592: **********************************************************************
593: SUBROUTINE BTPEC(N,PP,ISEED,JX)
594: BINOMIAL RANDOM VARIATE GENERATOR
595: MEAN .LT. 30 -- INVERSE CDF
596: MEAN .GE. 30 -- ALGORITHM BTPE: ACCEPTANCE-REJECTION VIA
597: FOUR REGION COMPOSITION. THE FOUR REGIONS ARE A TRIANGLE
598: (SYMMETRIC IN THE CENTER), A PAIR OF PARALLELOGRAMS (ABOVE
599: THE TRIANGLE), AND EXPONENTIAL LEFT AND RIGHT TAILS.
600: BTPE REFERS TO BINOMIAL-TRIANGLE-PARALLELOGRAM-EXPONENTIAL.
601: BTPEC REFERS TO BTPE AND "COMBINED." THUS BTPE IS THE
602: RESEARCH AND BTPEC IS THE IMPLEMENTATION OF A COMPLETE
603: USABLE ALGORITHM.
604: REFERENCE: VORATAS KACHITVICHYANUKUL AND BRUCE SCHMEISER,
605: "BINOMIAL RANDOM VARIATE GENERATION,"
606: COMMUNICATIONS OF THE ACM, FORTHCOMING
607: WRITTEN: SEPTEMBER 1980.
608: LAST REVISED: MAY 1985, JULY 1987
609: REQUIRED SUBPROGRAM: RAND() -- A UNIFORM (0,1) RANDOM NUMBER
610: GENERATOR
611: ARGUMENTS
612: N : NUMBER OF BERNOULLI TRIALS (INPUT)
613: PP : PROBABILITY OF SUCCESS IN EACH TRIAL (INPUT)
614: ISEED: RANDOM NUMBER SEED (INPUT AND OUTPUT)
615: JX: RANDOMLY GENERATED OBSERVATION (OUTPUT)
616: VARIABLES
617: PSAVE: VALUE OF PP FROM THE LAST CALL TO BTPEC
618: NSAVE: VALUE OF N FROM THE LAST CALL TO BTPEC
619: XNP: VALUE OF THE MEAN FROM THE LAST CALL TO BTPEC
620: P: PROBABILITY USED IN THE GENERATION PHASE OF BTPEC
621: FFM: TEMPORARY VARIABLE EQUAL TO XNP + P
622: M: INTEGER VALUE OF THE CURRENT MODE
623: FM: FLOATING POINT VALUE OF THE CURRENT MODE
624: XNPQ: TEMPORARY VARIABLE USED IN SETUP AND SQUEEZING STEPS
625: P1: AREA OF THE TRIANGLE
626: C: HEIGHT OF THE PARALLELOGRAMS
627: XM: CENTER OF THE TRIANGLE
628: XL: LEFT END OF THE TRIANGLE
629: XR: RIGHT END OF THE TRIANGLE
630: AL: TEMPORARY VARIABLE
631: XLL: RATE FOR THE LEFT EXPONENTIAL TAIL
632: XLR: RATE FOR THE RIGHT EXPONENTIAL TAIL
633: P2: AREA OF THE PARALLELOGRAMS
634: P3: AREA OF THE LEFT EXPONENTIAL TAIL
635: P4: AREA OF THE RIGHT EXPONENTIAL TAIL
636: U: A U(0,P4) RANDOM VARIATE USED FIRST TO SELECT ONE OF THE
637: FOUR REGIONS AND THEN CONDITIONALLY TO GENERATE A VALUE
638: FROM THE REGION
639: V: A U(0,1) RANDOM NUMBER USED TO GENERATE THE RANDOM VALUE
640: (REGION 1) OR TRANSFORMED INTO THE VARIATE TO ACCEPT OR
641: REJECT THE CANDIDATE VALUE
642: IX: INTEGER CANDIDATE VALUE
643: X: PRELIMINARY CONTINUOUS CANDIDATE VALUE IN REGION 2 LOGIC
644: AND A FLOATING POINT IX IN THE ACCEPT/REJECT LOGIC
645: K: ABSOLUTE VALUE OF (IX-M)
646: F: THE HEIGHT OF THE SCALED DENSITY FUNCTION USED IN THE
647: ACCEPT/REJECT DECISION WHEN BOTH M AND IX ARE SMALL
648: ALSO USED IN THE INVERSE TRANSFORMATION
649: R: THE RATIO P/Q
650: G: CONSTANT USED IN CALCULATION OF PROBABILITY
651: MP: MODE PLUS ONE, THE LOWER INDEX FOR EXPLICIT CALCULATION
652: OF F WHEN IX IS GREATER THAN M
653: IX1: CANDIDATE VALUE PLUS ONE, THE LOWER INDEX FOR EXPLICIT
654: CALCULATION OF F WHEN IX IS LESS THAN M
655: I: INDEX FOR EXPLICIT CALCULATION OF F FOR BTPE
656: AMAXP: MAXIMUM ERROR OF THE LOGARITHM OF NORMAL BOUND
657: YNORM: LOGARITHM OF NORMAL BOUND
658: ALV: NATURAL LOGARITHM OF THE ACCEPT/REJECT VARIATE V
659: X1,F1,Z,W,Z2,X2,F2, AND W2 ARE TEMPORARY VARIABLES TO BE
660: USED IN THE FINAL ACCEPT/REJECT TEST
661: QN: PROBABILITY OF NO SUCCESS IN N TRIALS
662: REMARK
663: IX AND JX COULD LOGICALLY BE THE SAME VARIABLE, WHICH WOULD
664: SAVE A MEMORY POSITION AND A LINE OF CODE. HOWEVER, SOME
665: COMPILERS (E.G.,CDC MNF) OPTIMIZE BETTER WHEN THE ARGUMENTS
666: ARE NOT INVOLVED.
667: ISEED NEEDS TO BE DOUBLE PRECISION IF THE IMSL ROUTINE
668: GGUBFS IS USED TO GENERATE UNIFORM RANDOM NUMBER, OTHERWISE
669: TYPE OF ISEED SHOULD BE DICTATED BY THE UNIFORM GENERATOR
670: **********************************************************************
671: *****DETERMINE APPROPRIATE ALGORITHM AND WHETHER SETUP IS NECESSARY
672: */
673: {
674: static float psave = -1.0;
675: static long nsave = -1;
676: static long ignbin,i,ix,ix1,k,m,mp,T1;
677: static float al,alv,amaxp,c,f,f1,f2,ffm,fm,g,p,p1,p2,p3,p4,q,qn,r,u,v,w,w2,x,x1,
678: x2,xl,xll,xlr,xm,xnp,xnpq,xr,ynorm,z,z2;
679:
680: if(pp != psave) goto S10;
681: if(n != nsave) goto S20;
682: if(xnp < 30.0) goto S150;
683: goto S30;
684: S10:
685: /*
686: *****SETUP, PERFORM ONLY WHEN PARAMETERS CHANGE
687: */
688: psave = pp;
689: p = min(psave,1.0-psave);
690: q = 1.0-p;
691: S20:
692: xnp = n*p;
693: nsave = n;
694: if(xnp < 30.0) goto S140;
695: ffm = xnp+p;
696: m = ffm;
697: fm = m;
698: xnpq = xnp*q;
699: p1 = (long) (2.195*sqrt(xnpq)-4.6*q)+0.5;
700: xm = fm+0.5;
701: xl = xm-p1;
702: xr = xm+p1;
703: c = 0.134+20.5/(15.3+fm);
704: al = (ffm-xl)/(ffm-xl*p);
705: xll = al*(1.0+0.5*al);
706: al = (xr-ffm)/(xr*q);
707: xlr = al*(1.0+0.5*al);
708: p2 = p1*(1.0+c+c);
709: p3 = p2+c/xll;
710: p4 = p3+c/xlr;
711: S30:
712: /*
713: *****GENERATE VARIATE
714: */
715: u = ranf()*p4;
716: v = ranf();
717: /*
718: TRIANGULAR REGION
719: */
720: if(u > p1) goto S40;
721: ix = xm-p1*v+u;
722: goto S170;
723: S40:
724: /*
725: PARALLELOGRAM REGION
726: */
727: if(u > p2) goto S50;
728: x = xl+(u-p1)/c;
729: v = v*c+1.0-abs(xm-x)/p1;
730: if(v > 1.0 || v <= 0.0) goto S30;
731: ix = x;
732: goto S70;
733: S50:
734: /*
735: LEFT TAIL
736: */
737: if(u > p3) goto S60;
738: ix = xl+log(v)/xll;
739: if(ix < 0) goto S30;
740: v *= ((u-p2)*xll);
741: goto S70;
742: S60:
743: /*
744: RIGHT TAIL
745: */
746: ix = xr-log(v)/xlr;
747: if(ix > n) goto S30;
748: v *= ((u-p3)*xlr);
749: S70:
750: /*
751: *****DETERMINE APPROPRIATE WAY TO PERFORM ACCEPT/REJECT TEST
752: */
753: k = abs(ix-m);
754: if(k > 20 && k < xnpq/2-1) goto S130;
755: /*
756: EXPLICIT EVALUATION
757: */
758: f = 1.0;
759: r = p/q;
760: g = (n+1)*r;
761: T1 = m-ix;
762: if(T1 < 0) goto S80;
763: else if(T1 == 0) goto S120;
764: else goto S100;
765: S80:
766: mp = m+1;
767: for(i=mp; i<=ix; i++) f *= (g/i-r);
768: goto S120;
769: S100:
770: ix1 = ix+1;
771: for(i=ix1; i<=m; i++) f /= (g/i-r);
772: S120:
773: if(v <= f) goto S170;
774: goto S30;
775: S130:
776: /*
777: SQUEEZING USING UPPER AND LOWER BOUNDS ON ALOG(F(X))
778: */
779: amaxp = k/xnpq*((k*(k/3.0+0.625)+0.1666666666666)/xnpq+0.5);
780: ynorm = -(k*k/(2.0*xnpq));
781: alv = log(v);
782: if(alv < ynorm-amaxp) goto S170;
783: if(alv > ynorm+amaxp) goto S30;
784: /*
785: STIRLING'S FORMULA TO MACHINE ACCURACY FOR
786: THE FINAL ACCEPTANCE/REJECTION TEST
787: */
788: x1 = ix+1.0;
789: f1 = fm+1.0;
790: z = n+1.0-fm;
791: w = n-ix+1.0;
792: z2 = z*z;
793: x2 = x1*x1;
794: f2 = f1*f1;
795: w2 = w*w;
796: if(alv <= xm*log(f1/x1)+(n-m+0.5)*log(z/w)+(ix-m)*log(w*p/(x1*q))+(13860.0-
797: (462.0-(132.0-(99.0-140.0/f2)/f2)/f2)/f2)/f1/166320.0+(13860.0-(462.0-
798: (132.0-(99.0-140.0/z2)/z2)/z2)/z2)/z/166320.0+(13860.0-(462.0-(132.0-
799: (99.0-140.0/x2)/x2)/x2)/x2)/x1/166320.0+(13860.0-(462.0-(132.0-(99.0
800: -140.0/w2)/w2)/w2)/w2)/w/166320.0) goto S170;
801: goto S30;
802: S140:
803: /*
804: INVERSE CDF LOGIC FOR MEAN LESS THAN 30
805: */
806: qn = pow(q,(double)n);
807: r = p/q;
808: g = r*(n+1);
809: S150:
810: ix = 0;
811: f = qn;
812: u = ranf();
813: S160:
814: if(u < f) goto S170;
815: if(ix > 110) goto S150;
816: u -= f;
817: ix += 1;
818: f *= (g/ix-r);
819: goto S160;
820: S170:
821: if(psave > 0.5) ix = n-ix;
822: ignbin = ix;
823: return ignbin;
824: }
825: long ignpoi(float mu)
826: /*
827: **********************************************************************
828: long ignpoi(float mu)
829: GENerate POIsson random deviate
830: Function
831: Generates a single random deviate from a Poisson
832: distribution with mean AV.
833: Arguments
834: av --> The mean of the Poisson distribution from which
835: a random deviate is to be generated.
836: genexp <-- The random deviate.
837: Method
838: Renames KPOIS from TOMS as slightly modified by BWB to use RANF
839: instead of SUNIF.
840: For details see:
841: Ahrens, J.H. and Dieter, U.
842: Computer Generation of Poisson Deviates
843: From Modified Normal Distributions.
844: ACM Trans. Math. Software, 8, 2
845: (June 1982),163-179
846: **********************************************************************
847: **********************************************************************
848:
849:
850: P O I S S O N DISTRIBUTION
851:
852:
853: **********************************************************************
854: **********************************************************************
855:
856: FOR DETAILS SEE:
857:
858: AHRENS, J.H. AND DIETER, U.
859: COMPUTER GENERATION OF POISSON DEVIATES
860: FROM MODIFIED NORMAL DISTRIBUTIONS.
861: ACM TRANS. MATH. SOFTWARE, 8,2 (JUNE 1982), 163 - 179.
862:
863: (SLIGHTLY MODIFIED VERSION OF THE PROGRAM IN THE ABOVE ARTICLE)
864:
865: **********************************************************************
866: INTEGER FUNCTION IGNPOI(IR,MU)
867: INPUT: IR=CURRENT STATE OF BASIC RANDOM NUMBER GENERATOR
868: MU=MEAN MU OF THE POISSON DISTRIBUTION
869: OUTPUT: IGNPOI=SAMPLE FROM THE POISSON-(MU)-DISTRIBUTION
870: MUPREV=PREVIOUS MU, MUOLD=MU AT LAST EXECUTION OF STEP P OR B.
871: TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT
872: COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL
873: SEPARATION OF CASES A AND B
874: */
875: {
876: extern float fsign( float num, float sign );
877: static float a0 = -0.5;
878: static float a1 = 0.3333333;
879: static float a2 = -0.2500068;
880: static float a3 = 0.2000118;
881: static float a4 = -0.1661269;
882: static float a5 = 0.1421878;
883: static float a6 = -0.1384794;
884: static float a7 = 0.125006;
885: static float muold = 0.0;
886: static float muprev = 0.0;
887: static float fact[10] = {
888: 1.0,1.0,2.0,6.0,24.0,120.0,720.0,5040.0,40320.0,362880.0
889: };
890: static long ignpoi,j,k,kflag,l,m;
891: static float b1,b2,c,c0,c1,c2,c3,d,del,difmuk,e,fk,fx,fy,g,omega,p,p0,px,py,q,s,
892: t,u,v,x,xx,pp[35];
893:
894: if(mu == muprev) goto S10;
895: if(mu < 10.0) goto S120;
896: /*
897: C A S E A. (RECALCULATION OF S,D,L IF MU HAS CHANGED)
898: */
899: muprev = mu;
900: s = sqrt(mu);
901: d = 6.0*mu*mu;
902: /*
903: THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL
904: PROBABILITIES FK WHENEVER K >= M(MU). L=IFIX(MU-1.1484)
905: IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 .
906: */
907: l = (long) (mu-1.1484);
908: S10:
909: /*
910: STEP N. NORMAL SAMPLE - SNORM(IR) FOR STANDARD NORMAL DEVIATE
911: */
912: g = mu+s*snorm();
913: if(g < 0.0) goto S20;
914: ignpoi = (long) (g);
915: /*
916: STEP I. IMMEDIATE ACCEPTANCE IF IGNPOI IS LARGE ENOUGH
917: */
918: if(ignpoi >= l) return ignpoi;
919: /*
920: STEP S. SQUEEZE ACCEPTANCE - SUNIF(IR) FOR (0,1)-SAMPLE U
921: */
922: fk = (float)ignpoi;
923: difmuk = mu-fk;
924: u = ranf();
925: if(d*u >= difmuk*difmuk*difmuk) return ignpoi;
926: S20:
927: /*
928: STEP P. PREPARATIONS FOR STEPS Q AND H.
929: (RECALCULATIONS OF PARAMETERS IF NECESSARY)
930: .3989423=(2*PI)**(-.5) .416667E-1=1./24. .1428571=1./7.
931: THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE
932: APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK.
933: C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION.
934: */
935: if(mu == muold) goto S30;
936: muold = mu;
937: omega = 0.3989423/s;
938: b1 = 4.166667E-2/mu;
939: b2 = 0.3*b1*b1;
940: c3 = 0.1428571*b1*b2;
941: c2 = b2-15.0*c3;
942: c1 = b1-6.0*b2+45.0*c3;
943: c0 = 1.0-b1+3.0*b2-15.0*c3;
944: c = 0.1069/mu;
945: S30:
946: if(g < 0.0) goto S50;
947: /*
948: 'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN)
949: */
950: kflag = 0;
951: goto S70;
952: S40:
953: /*
954: STEP Q. QUOTIENT ACCEPTANCE (RARE CASE)
955: */
956: if(fy-u*fy <= py*exp(px-fx)) return ignpoi;
957: S50:
958: /*
959: STEP E. EXPONENTIAL SAMPLE - SEXPO(IR) FOR STANDARD EXPONENTIAL
960: DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT'
961: (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.)
962: */
963: e = sexpo();
964: u = ranf();
965: u += (u-1.0);
966: t = 1.8+fsign(e,u);
967: if(t <= -0.6744) goto S50;
968: ignpoi = (long) (mu+s*t);
969: fk = (float)ignpoi;
970: difmuk = mu-fk;
971: /*
972: 'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN)
973: */
974: kflag = 1;
975: goto S70;
976: S60:
977: /*
978: STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION)
979: */
980: if(c*fabs(u) > py*exp(px+e)-fy*exp(fx+e)) goto S50;
981: return ignpoi;
982: S70:
983: /*
984: STEP F. 'SUBROUTINE' F. CALCULATION OF PX,PY,FX,FY.
985: CASE IGNPOI .LT. 10 USES FACTORIALS FROM TABLE FACT
986: */
987: if(ignpoi >= 10) goto S80;
988: px = -mu;
989: py = pow(mu,(double)ignpoi)/ *(fact+ignpoi);
990: goto S110;
991: S80:
992: /*
993: CASE IGNPOI .GE. 10 USES POLYNOMIAL APPROXIMATION
994: A0-A7 FOR ACCURACY WHEN ADVISABLE
995: .8333333E-1=1./12. .3989423=(2*PI)**(-.5)
996: */
997: del = 8.333333E-2/fk;
998: del -= (4.8*del*del*del);
999: v = difmuk/fk;
1000: if(fabs(v) <= 0.25) goto S90;
1001: px = fk*log(1.0+v)-difmuk-del;
1002: goto S100;
1003: S90:
1004: px = fk*v*v*(((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v+a0)-del;
1005: S100:
1006: py = 0.3989423/sqrt(fk);
1007: S110:
1008: x = (0.5-difmuk)/s;
1009: xx = x*x;
1010: fx = -0.5*xx;
1011: fy = omega*(((c3*xx+c2)*xx+c1)*xx+c0);
1012: if(kflag <= 0) goto S40;
1013: goto S60;
1014: S120:
1015: /*
1016: C A S E B. (START NEW TABLE AND CALCULATE P0 IF NECESSARY)
1017: */
1018: muprev = 0.0;
1019: if(mu == muold) goto S130;
1020: muold = mu;
1021: m = max(1L,(long) (mu));
1022: l = 0;
1023: p = exp(-mu);
1024: q = p0 = p;
1025: S130:
1026: /*
1027: STEP U. UNIFORM SAMPLE FOR INVERSION METHOD
1028: */
1029: u = ranf();
1030: ignpoi = 0;
1031: if(u <= p0) return ignpoi;
1032: /*
1033: STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE
1034: PP-TABLE OF CUMULATIVE POISSON PROBABILITIES
1035: (0.458=PP(9) FOR MU=10)
1036: */
1037: if(l == 0) goto S150;
1038: j = 1;
1039: if(u > 0.458) j = min(l,m);
1040: for(k=j; k<=l; k++) {
1041: if(u <= *(pp+k-1)) goto S180;
1042: }
1043: if(l == 35) goto S130;
1044: S150:
1045: /*
1046: STEP C. CREATION OF NEW POISSON PROBABILITIES P
1047: AND THEIR CUMULATIVES Q=PP(K)
1048: */
1049: l += 1;
1050: for(k=l; k<=35; k++) {
1051: p = p*mu/(float)k;
1052: q += p;
1053: *(pp+k-1) = q;
1054: if(u <= q) goto S170;
1055: }
1056: l = 35;
1057: goto S130;
1058: S170:
1059: l = k;
1060: S180:
1061: ignpoi = k;
1062: return ignpoi;
1063: }
1064: long ignuin(long low,long high)
1065: /*
1066: **********************************************************************
1067: long ignuin(long low,long high)
1068: GeNerate Uniform INteger
1069: Function
1070: Generates an integer uniformly distributed between LOW and HIGH.
1071: Arguments
1072: low --> Low bound (inclusive) on integer value to be generated
1073: high --> High bound (inclusive) on integer value to be generated
1074: Note
1075: If (HIGH-LOW) > 2,147,483,561 prints error message on * unit and
1076: stops the program.
1077: **********************************************************************
1078: IGNLGI generates integers between 1 and 2147483562
1079: MAXNUM is 1 less than maximum generable value
1080: */
1081: {
1082: #define maxnum 2147483561L
1083: static long ignuin,ign,maxnow,range,ranp1;
1084:
1085: if(!(low > high)) goto S10;
1086: puts(" low > high in ignuin - ABORT");
1087: exit(1);
1088:
1089: S10:
1090: range = high-low;
1091: if(!(range > maxnum)) goto S20;
1092: puts(" high - low too large in ignuin - ABORT");
1093: exit(1);
1094:
1095: S20:
1096: if(!(low == high)) goto S30;
1097: ignuin = low;
1098: return ignuin;
1099:
1100: S30:
1101: /*
1102: Number to be generated should be in range 0..RANGE
1103: Set MAXNOW so that the number of integers in 0..MAXNOW is an
1104: integral multiple of the number in 0..RANGE
1105: */
1106: ranp1 = range+1;
1107: maxnow = maxnum/ranp1*ranp1;
1108: S40:
1109: ign = ignlgi()-1;
1110: if(!(ign <= maxnow)) goto S50;
1111: ignuin = low+ign%ranp1;
1112: return ignuin;
1113: S50:
1114: goto S40;
1115: #undef maxnum
1116: #undef err1
1117: #undef err2
1118: }
1119: long lennob( char *str )
1120: /*
1121: Returns the length of str ignoring trailing blanks but not
1122: other white space.
1123: */
1124: {
1125: long i, i_nb;
1126:
1127: for (i=0, i_nb= -1L; *(str+i); i++)
1128: if ( *(str+i) != ' ' ) i_nb = i;
1129: return (i_nb+1);
1130: }
1131: long mltmod(long a,long s,long m)
1132: /*
1133: **********************************************************************
1134: long mltmod(long a,long s,long m)
1135: Returns (A*S) MOD M
1136: This is a transcription from Pascal to Fortran of routine
1137: MULtMod_Decompos from the paper
1138: L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
1139: with Splitting Facilities." ACM Transactions on Mathematical
1140: Software, 17:98-111 (1991)
1141: Arguments
1142: a, s, m -->
1143: **********************************************************************
1144: */
1145: {
1146: #define h 32768L
1147: static long mltmod,a0,a1,k,p,q,qh,rh;
1148: /*
1149: H = 2**((b-2)/2) where b = 32 because we are using a 32 bit
1150: machine. On a different machine recompute H
1151: */
1152: if(!(a <= 0 || a >= m || s <= 0 || s >= m)) goto S10;
1153: puts(" a, m, s out of order in mltmod - ABORT!");
1154: printf(" a = %12ld s = %12ld m = %12ld\n",a,s,m);
1155: puts(" mltmod requires: 0 < a < m; 0 < s < m");
1156: exit(1);
1157: S10:
1158: if(!(a < h)) goto S20;
1159: a0 = a;
1160: p = 0;
1161: goto S120;
1162: S20:
1163: a1 = a/h;
1164: a0 = a-h*a1;
1165: qh = m/h;
1166: rh = m-h*qh;
1167: if(!(a1 >= h)) goto S50;
1168: a1 -= h;
1169: k = s/qh;
1170: p = h*(s-k*qh)-k*rh;
1171: S30:
1172: if(!(p < 0)) goto S40;
1173: p += m;
1174: goto S30;
1175: S40:
1176: goto S60;
1177: S50:
1178: p = 0;
1179: S60:
1180: /*
1181: P = (A2*S*H)MOD M
1182: */
1183: if(!(a1 != 0)) goto S90;
1184: q = m/a1;
1185: k = s/q;
1186: p -= (k*(m-a1*q));
1187: if(p > 0) p -= m;
1188: p += (a1*(s-k*q));
1189: S70:
1190: if(!(p < 0)) goto S80;
1191: p += m;
1192: goto S70;
1193: S90:
1194: S80:
1195: k = p/qh;
1196: /*
1197: P = ((A2*H + A1)*S)MOD M
1198: */
1199: p = h*(p-k*qh)-k*rh;
1200: S100:
1201: if(!(p < 0)) goto S110;
1202: p += m;
1203: goto S100;
1204: S120:
1205: S110:
1206: if(!(a0 != 0)) goto S150;
1207: /*
1208: P = ((A2*H + A1)*H*S)MOD M
1209: */
1210: q = m/a0;
1211: k = s/q;
1212: p -= (k*(m-a0*q));
1213: if(p > 0) p -= m;
1214: p += (a0*(s-k*q));
1215: S130:
1216: if(!(p < 0)) goto S140;
1217: p += m;
1218: goto S130;
1219: S150:
1220: S140:
1221: mltmod = p;
1222: return mltmod;
1223: #undef h
1224: }
1225: void phrtsd(char* phrase,long *seed1,long *seed2)
1226: /*
1227: **********************************************************************
1228: void phrtsd(char* phrase,long *seed1,long *seed2)
1229: PHRase To SeeDs
1230:
1231: Function
1232:
1233: Uses a phrase (character string) to generate two seeds for the RGN
1234: random number generator.
1235: Arguments
1236: phrase --> Phrase to be used for random number generation
1237:
1238: seed1 <-- First seed for generator
1239:
1240: seed2 <-- Second seed for generator
1241:
1242: Note
1243:
1244: Trailing blanks are eliminated before the seeds are generated.
1245: Generated seed values will fall in the range 1..2^30
1246: (1..1,073,741,824)
1247: **********************************************************************
1248: */
1249: {
1250:
1251: static char table[] =
1252: "abcdefghijklmnopqrstuvwxyz\
1253: ABCDEFGHIJKLMNOPQRSTUVWXYZ\
1254: 0123456789\
1255: !@#$%^&*()_+[];:'\\\"<>?,./";
1256:
1257: long ix;
1258:
1259: static long twop30 = 1073741824L;
1260: static long shift[5] = {
1261: 1L,64L,4096L,262144L,16777216L
1262: };
1263: static long i,ichr,j,lphr,values[5];
1264: extern long lennob(char *str);
1265:
1266: *seed1 = 1234567890L;
1267: *seed2 = 123456789L;
1268: lphr = lennob(phrase);
1269: if(lphr < 1) return;
1270: for(i=0; i<=(lphr-1); i++) {
1271: for (ix=0; table[ix]; ix++) if (*(phrase+i) == table[ix]) break;
1272: if (!table[ix]) ix = 0;
1273: ichr = ix % 64;
1274: if(ichr == 0) ichr = 63;
1275: for(j=1; j<=5; j++) {
1276: *(values+j-1) = ichr-j;
1277: if(*(values+j-1) < 1) *(values+j-1) += 63;
1278: }
1279: for(j=1; j<=5; j++) {
1280: *seed1 = ( *seed1+*(shift+j-1)**(values+j-1) ) % twop30;
1281: *seed2 = ( *seed2+*(shift+j-1)**(values+6-j-1) ) % twop30;
1282: }
1283: }
1284: #undef twop30
1285: }
1286: float ranf(void)
1287: /*
1288: **********************************************************************
1289: float ranf(void)
1290: RANDom number generator as a Function
1291: Returns a random floating point number from a uniform distribution
1292: over 0 - 1 (endpoints of this interval are not returned) using the
1293: current generator
1294: This is a transcription from Pascal to Fortran of routine
1295: Uniform_01 from the paper
1296: L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
1297: with Splitting Facilities." ACM Transactions on Mathematical
1298: Software, 17:98-111 (1991)
1299: **********************************************************************
1300: */
1301: {
1302: static float ranf;
1303: /*
1304: 4.656613057E-10 is 1/M1 M1 is set in a data statement in IGNLGI
1305: and is currently 2147483563. If M1 changes, change this also.
1306: */
1307: ranf = ignlgi()*4.656613057E-10;
1308: return ranf;
1309: }
1310: void setgmn(float *meanv,float *covm,long p,float *parm)
1311: /*
1312: **********************************************************************
1313: void setgmn(float *meanv,float *covm,long p,float *parm)
1314: SET Generate Multivariate Normal random deviate
1315: Function
1316: Places P, MEANV, and the Cholesky factoriztion of COVM
1317: in GENMN.
1318: Arguments
1319: meanv --> Mean vector of multivariate normal distribution.
1320: covm <--> (Input) Covariance matrix of the multivariate
1321: normal distribution
1322: (Output) Destroyed on output
1323: p --> Dimension of the normal, or length of MEANV.
1324: parm <-- Array of parameters needed to generate multivariate norma
1325: deviates (P, MEANV and Cholesky decomposition of
1326: COVM).
1327: 1 : 1 - P
1328: 2 : P + 1 - MEANV
1329: P+2 : P*(P+3)/2 + 1 - Cholesky decomposition of COVM
1330: Needed dimension is (p*(p+3)/2 + 1)
1331: **********************************************************************
1332: */
1333: {
1334: extern void spofa(float *a,long lda,long n,long *info);
1335: static long T1;
1336: static long i,icount,info,j,D2,D3,D4,D5;
1337: T1 = p*(p+3)/2+1;
1338: /*
1339: TEST THE INPUT
1340: */
1341: if(!(p <= 0)) goto S10;
1342: puts("P nonpositive in SETGMN");
1343: printf("Value of P: %12ld\n",p);
1344: exit(1);
1345: S10:
1346: *parm = p;
1347: /*
1348: PUT P AND MEANV INTO PARM
1349: */
1350: for(i=2,D2=1,D3=(p+1-i+D2)/D2; D3>0; D3--,i+=D2) *(parm+i-1) = *(meanv+i-2);
1351: /*
1352: Cholesky decomposition to find A s.t. trans(A)*(A) = COVM
1353: */
1354: spofa(covm,p,p,&info);
1355: if(!(info != 0)) goto S30;
1356: puts(" COVM not positive definite in SETGMN");
1357: exit(1);
1358: S30:
1359: icount = p+1;
1360: /*
1361: PUT UPPER HALF OF A, WHICH IS NOW THE CHOLESKY FACTOR, INTO PARM
1362: COVM(1,1) = PARM(P+2)
1363: COVM(1,2) = PARM(P+3)
1364: :
1365: COVM(1,P) = PARM(2P+1)
1366: COVM(2,2) = PARM(2P+2) ...
1367: */
1368: for(i=1,D4=1,D5=(p-i+D4)/D4; D5>0; D5--,i+=D4) {
1369: for(j=i-1; j<p; j++) {
1370: icount += 1;
1371: *(parm+icount-1) = *(covm+i-1+j*p);
1372: }
1373: }
1374: }
1375: float sexpo(void)
1376: /*
1377: **********************************************************************
1378:
1379:
1380: (STANDARD-) E X P O N E N T I A L DISTRIBUTION
1381:
1382:
1383: **********************************************************************
1384: **********************************************************************
1385:
1386: FOR DETAILS SEE:
1387:
1388: AHRENS, J.H. AND DIETER, U.
1389: COMPUTER METHODS FOR SAMPLING FROM THE
1390: EXPONENTIAL AND NORMAL DISTRIBUTIONS.
1391: COMM. ACM, 15,10 (OCT. 1972), 873 - 882.
1392:
1393: ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM
1394: 'SA' IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION)
1395:
1396: Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of
1397: SUNIF. The argument IR thus goes away.
1398:
1399: **********************************************************************
1400: Q(N) = SUM(ALOG(2.0)**K/K!) K=1,..,N , THE HIGHEST N
1401: (HERE 8) IS DETERMINED BY Q(N)=1.0 WITHIN STANDARD PRECISION
1402: */
1403: {
1404: static float q[8] = {
1405: 0.6931472,0.9333737,0.9888778,0.9984959,0.9998293,0.9999833,0.9999986,1.0
1406: };
1407: static long i;
1408: static float sexpo,a,u,ustar,umin;
1409: static float *q1 = q;
1410: a = 0.0;
1411: u = ranf();
1412: goto S30;
1413: S20:
1414: a += *q1;
1415: S30:
1416: u += u;
1417: if(u <= 1.0) goto S20;
1418: u -= 1.0;
1419: if(u > *q1) goto S60;
1420: sexpo = a+u;
1421: return sexpo;
1422: S60:
1423: i = 1;
1424: ustar = ranf();
1425: umin = ustar;
1426: S70:
1427: ustar = ranf();
1428: if(ustar < umin) umin = ustar;
1429: i += 1;
1430: if(u > *(q+i-1)) goto S70;
1431: sexpo = a+umin**q1;
1432: return sexpo;
1433: }
1434: float sgamma(float a)
1435: /*
1436: **********************************************************************
1437:
1438:
1439: (STANDARD-) G A M M A DISTRIBUTION
1440:
1441:
1442: **********************************************************************
1443: **********************************************************************
1444:
1445: PARAMETER A >= 1.0 !
1446:
1447: **********************************************************************
1448:
1449: FOR DETAILS SEE:
1450:
1451: AHRENS, J.H. AND DIETER, U.
1452: GENERATING GAMMA VARIATES BY A
1453: MODIFIED REJECTION TECHNIQUE.
1454: COMM. ACM, 25,1 (JAN. 1982), 47 - 54.
1455:
1456: STEP NUMBERS CORRESPOND TO ALGORITHM 'GD' IN THE ABOVE PAPER
1457: (STRAIGHTFORWARD IMPLEMENTATION)
1458:
1459: Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of
1460: SUNIF. The argument IR thus goes away.
1461:
1462: **********************************************************************
1463:
1464: PARAMETER 0.0 < A < 1.0 !
1465:
1466: **********************************************************************
1467:
1468: FOR DETAILS SEE:
1469:
1470: AHRENS, J.H. AND DIETER, U.
1471: COMPUTER METHODS FOR SAMPLING FROM GAMMA,
1472: BETA, POISSON AND BINOMIAL DISTRIBUTIONS.
1473: COMPUTING, 12 (1974), 223 - 246.
1474:
1475: (ADAPTED IMPLEMENTATION OF ALGORITHM 'GS' IN THE ABOVE PAPER)
1476:
1477: **********************************************************************
1478: INPUT: A =PARAMETER (MEAN) OF THE STANDARD GAMMA DISTRIBUTION
1479: OUTPUT: SGAMMA = SAMPLE FROM THE GAMMA-(A)-DISTRIBUTION
1480: COEFFICIENTS Q(K) - FOR Q0 = SUM(Q(K)*A**(-K))
1481: COEFFICIENTS A(K) - FOR Q = Q0+(T*T/2)*SUM(A(K)*V**K)
1482: COEFFICIENTS E(K) - FOR EXP(Q)-1 = SUM(E(K)*Q**K)
1483: PREVIOUS A PRE-SET TO ZERO - AA IS A', AAA IS A"
1484: SQRT32 IS THE SQUAREROOT OF 32 = 5.656854249492380
1485: */
1486: {
1487: extern float fsign( float num, float sign );
1488: static float q1 = 4.166669E-2;
1489: static float q2 = 2.083148E-2;
1490: static float q3 = 8.01191E-3;
1491: static float q4 = 1.44121E-3;
1492: static float q5 = -7.388E-5;
1493: static float q6 = 2.4511E-4;
1494: static float q7 = 2.424E-4;
1495: static float a1 = 0.3333333;
1496: static float a2 = -0.250003;
1497: static float a3 = 0.2000062;
1498: static float a4 = -0.1662921;
1499: static float a5 = 0.1423657;
1500: static float a6 = -0.1367177;
1501: static float a7 = 0.1233795;
1502: static float e1 = 1.0;
1503: static float e2 = 0.4999897;
1504: static float e3 = 0.166829;
1505: static float e4 = 4.07753E-2;
1506: static float e5 = 1.0293E-2;
1507: static float aa = 0.0;
1508: static float aaa = 0.0;
1509: static float sqrt32 = 5.656854;
1510: static float sgamma,s2,s,d,t,x,u,r,q0,b,si,c,v,q,e,w,p;
1511: if(a == aa) goto S10;
1512: if(a < 1.0) goto S120;
1513: /*
1514: STEP 1: RECALCULATIONS OF S2,S,D IF A HAS CHANGED
1515: */
1516: aa = a;
1517: s2 = a-0.5;
1518: s = sqrt(s2);
1519: d = sqrt32-12.0*s;
1520: S10:
1521: /*
1522: STEP 2: T=STANDARD NORMAL DEVIATE,
1523: X=(S,1/2)-NORMAL DEVIATE.
1524: IMMEDIATE ACCEPTANCE (I)
1525: */
1526: t = snorm();
1527: x = s+0.5*t;
1528: sgamma = x*x;
1529: if(t >= 0.0) return sgamma;
1530: /*
1531: STEP 3: U= 0,1 -UNIFORM SAMPLE. SQUEEZE ACCEPTANCE (S)
1532: */
1533: u = ranf();
1534: if(d*u <= t*t*t) return sgamma;
1535: /*
1536: STEP 4: RECALCULATIONS OF Q0,B,SI,C IF NECESSARY
1537: */
1538: if(a == aaa) goto S40;
1539: aaa = a;
1540: r = 1.0/ a;
1541: q0 = ((((((q7*r+q6)*r+q5)*r+q4)*r+q3)*r+q2)*r+q1)*r;
1542: /*
1543: APPROXIMATION DEPENDING ON SIZE OF PARAMETER A
1544: THE CONSTANTS IN THE EXPRESSIONS FOR B, SI AND
1545: C WERE ESTABLISHED BY NUMERICAL EXPERIMENTS
1546: */
1547: if(a <= 3.686) goto S30;
1548: if(a <= 13.022) goto S20;
1549: /*
1550: CASE 3: A .GT. 13.022
1551: */
1552: b = 1.77;
1553: si = 0.75;
1554: c = 0.1515/s;
1555: goto S40;
1556: S20:
1557: /*
1558: CASE 2: 3.686 .LT. A .LE. 13.022
1559: */
1560: b = 1.654+7.6E-3*s2;
1561: si = 1.68/s+0.275;
1562: c = 6.2E-2/s+2.4E-2;
1563: goto S40;
1564: S30:
1565: /*
1566: CASE 1: A .LE. 3.686
1567: */
1568: b = 0.463+s-0.178*s2;
1569: si = 1.235;
1570: c = 0.195/s-7.9E-2+1.6E-2*s;
1571: S40:
1572: /*
1573: STEP 5: NO QUOTIENT TEST IF X NOT POSITIVE
1574: */
1575: if(x <= 0.0) goto S70;
1576: /*
1577: STEP 6: CALCULATION OF V AND QUOTIENT Q
1578: */
1579: v = t/(s+s);
1580: if(fabs(v) <= 0.25) goto S50;
1581: q = q0-s*t+0.25*t*t+(s2+s2)*log(1.0+v);
1582: goto S60;
1583: S50:
1584: q = q0+0.5*t*t*((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v;
1585: S60:
1586: /*
1587: STEP 7: QUOTIENT ACCEPTANCE (Q)
1588: */
1589: if(log(1.0-u) <= q) return sgamma;
1590: S70:
1591: /*
1592: STEP 8: E=STANDARD EXPONENTIAL DEVIATE
1593: U= 0,1 -UNIFORM DEVIATE
1594: T=(B,SI)-DOUBLE EXPONENTIAL (LAPLACE) SAMPLE
1595: */
1596: e = sexpo();
1597: u = ranf();
1598: u += (u-1.0);
1599: t = b+fsign(si*e,u);
1600: /*
1601: STEP 9: REJECTION IF T .LT. TAU(1) = -.71874483771719
1602: */
1603: if(t < -0.7187449) goto S70;
1604: /*
1605: STEP 10: CALCULATION OF V AND QUOTIENT Q
1606: */
1607: v = t/(s+s);
1608: if(fabs(v) <= 0.25) goto S80;
1609: q = q0-s*t+0.25*t*t+(s2+s2)*log(1.0+v);
1610: goto S90;
1611: S80:
1612: q = q0+0.5*t*t*((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v;
1613: S90:
1614: /*
1615: STEP 11: HAT ACCEPTANCE (H) (IF Q NOT POSITIVE GO TO STEP 8)
1616: */
1617: if(q <= 0.0) goto S70;
1618: if(q <= 0.5) goto S100;
1619: w = exp(q)-1.0;
1620: goto S110;
1621: S100:
1622: w = ((((e5*q+e4)*q+e3)*q+e2)*q+e1)*q;
1623: S110:
1624: /*
1625: IF T IS REJECTED, SAMPLE AGAIN AT STEP 8
1626: */
1627: if(c*fabs(u) > w*exp(e-0.5*t*t)) goto S70;
1628: x = s+0.5*t;
1629: sgamma = x*x;
1630: return sgamma;
1631: S120:
1632: /*
1633: ALTERNATE METHOD FOR PARAMETERS A BELOW 1 (.3678794=EXP(-1.))
1634: */
1635: aa = 0.0;
1636: b = 1.0+0.3678794*a;
1637: S130:
1638: p = b*ranf();
1639: if(p >= 1.0) goto S140;
1640: sgamma = exp(log(p)/ a);
1641: if(sexpo() < sgamma) goto S130;
1642: return sgamma;
1643: S140:
1644: sgamma = -log((b-p)/ a);
1645: if(sexpo() < (1.0-a)*log(sgamma)) goto S130;
1646: return sgamma;
1647: }
1648: float snorm(void)
1649: /*
1650: **********************************************************************
1651:
1652:
1653: (STANDARD-) N O R M A L DISTRIBUTION
1654:
1655:
1656: **********************************************************************
1657: **********************************************************************
1658:
1659: FOR DETAILS SEE:
1660:
1661: AHRENS, J.H. AND DIETER, U.
1662: EXTENSIONS OF FORSYTHE'S METHOD FOR RANDOM
1663: SAMPLING FROM THE NORMAL DISTRIBUTION.
1664: MATH. COMPUT., 27,124 (OCT. 1973), 927 - 937.
1665:
1666: ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM 'FL'
1667: (M=5) IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION)
1668:
1669: Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of
1670: SUNIF. The argument IR thus goes away.
1671:
1672: **********************************************************************
1673: THE DEFINITIONS OF THE CONSTANTS A(K), D(K), T(K) AND
1674: H(K) ARE ACCORDING TO THE ABOVEMENTIONED ARTICLE
1675: */
1676: {
1677: static float a[32] = {
1678: 0.0,3.917609E-2,7.841241E-2,0.11777,0.1573107,0.1970991,0.2372021,0.2776904,
1679: 0.3186394,0.36013,0.4022501,0.4450965,0.4887764,0.5334097,0.5791322,
1680: 0.626099,0.6744898,0.7245144,0.7764218,0.8305109,0.8871466,0.9467818,
1681: 1.00999,1.077516,1.150349,1.229859,1.318011,1.417797,1.534121,1.67594,
1682: 1.862732,2.153875
1683: };
1684: static float d[31] = {
1685: 0.0,0.0,0.0,0.0,0.0,0.2636843,0.2425085,0.2255674,0.2116342,0.1999243,
1686: 0.1899108,0.1812252,0.1736014,0.1668419,0.1607967,0.1553497,0.1504094,
1687: 0.1459026,0.14177,0.1379632,0.1344418,0.1311722,0.128126,0.1252791,
1688: 0.1226109,0.1201036,0.1177417,0.1155119,0.1134023,0.1114027,0.1095039
1689: };
1690: static float t[31] = {
1691: 7.673828E-4,2.30687E-3,3.860618E-3,5.438454E-3,7.0507E-3,8.708396E-3,
1692: 1.042357E-2,1.220953E-2,1.408125E-2,1.605579E-2,1.81529E-2,2.039573E-2,
1693: 2.281177E-2,2.543407E-2,2.830296E-2,3.146822E-2,3.499233E-2,3.895483E-2,
1694: 4.345878E-2,4.864035E-2,5.468334E-2,6.184222E-2,7.047983E-2,8.113195E-2,
1695: 9.462444E-2,0.1123001,0.136498,0.1716886,0.2276241,0.330498,0.5847031
1696: };
1697: static float h[31] = {
1698: 3.920617E-2,3.932705E-2,3.951E-2,3.975703E-2,4.007093E-2,4.045533E-2,
1699: 4.091481E-2,4.145507E-2,4.208311E-2,4.280748E-2,4.363863E-2,4.458932E-2,
1700: 4.567523E-2,4.691571E-2,4.833487E-2,4.996298E-2,5.183859E-2,5.401138E-2,
1701: 5.654656E-2,5.95313E-2,6.308489E-2,6.737503E-2,7.264544E-2,7.926471E-2,
1702: 8.781922E-2,9.930398E-2,0.11556,0.1404344,0.1836142,0.2790016,0.7010474
1703: };
1704: static long i;
1705: static float snorm,u,s,ustar,aa,w,y,tt;
1706: u = ranf();
1707: s = 0.0;
1708: if(u > 0.5) s = 1.0;
1709: u += (u-s);
1710: u = 32.0*u;
1711: i = (long) (u);
1712: if(i == 32) i = 31;
1713: if(i == 0) goto S100;
1714: /*
1715: START CENTER
1716: */
1717: ustar = u-(float)i;
1718: aa = *(a+i-1);
1719: S40:
1720: if(ustar <= *(t+i-1)) goto S60;
1721: w = (ustar-*(t+i-1))**(h+i-1);
1722: S50:
1723: /*
1724: EXIT (BOTH CASES)
1725: */
1726: y = aa+w;
1727: snorm = y;
1728: if(s == 1.0) snorm = -y;
1729: return snorm;
1730: S60:
1731: /*
1732: CENTER CONTINUED
1733: */
1734: u = ranf();
1735: w = u*(*(a+i)-aa);
1736: tt = (0.5*w+aa)*w;
1737: goto S80;
1738: S70:
1739: tt = u;
1740: ustar = ranf();
1741: S80:
1742: if(ustar > tt) goto S50;
1743: u = ranf();
1744: if(ustar >= u) goto S70;
1745: ustar = ranf();
1746: goto S40;
1747: S100:
1748: /*
1749: START TAIL
1750: */
1751: i = 6;
1752: aa = *(a+31);
1753: goto S120;
1754: S110:
1755: aa += *(d+i-1);
1756: i += 1;
1757: S120:
1758: u += u;
1759: if(u < 1.0) goto S110;
1760: u -= 1.0;
1761: S140:
1762: w = u**(d+i-1);
1763: tt = (0.5*w+aa)*w;
1764: goto S160;
1765: S150:
1766: tt = u;
1767: S160:
1768: ustar = ranf();
1769: if(ustar > tt) goto S50;
1770: u = ranf();
1771: if(ustar >= u) goto S150;
1772: u = ranf();
1773: goto S140;
1774: }
1775: float fsign( float num, float sign )
1776: /* Transfers sign of argument sign to argument num */
1777: {
1778: if ( ( sign>0.0f && num<0.0f ) || ( sign<0.0f && num>0.0f ) )
1779: return -num;
1780: else return num;
1781: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/back.gif){kind=icon}
Return to ranlib.c CVS log ![[TXT]](/icons/text.gif){kind=icon}
![[DIR]](/icons/dir.gif){kind=icon}
Up to [LON-CAPA] / capa / capa51 / pProj