1 SUBROUTINE stat(x,n,av,var,xmin,xmax)
35 IF (x(i).LT.xmin) xmin = x(i)
36 IF (x(i).GT.xmax) xmax = x(i)
41 sum = sum + (x(i)-av)**2
52 parameter(mxwh=15,mxncat=100)
55 REAL av,avtr,var,vartr,
xmin,
xmax,pevt,psum,rtry
56 INTEGER i,is1,is2,itmp,iwhich,j,mxint,nperm,nrep,ntot,ntry,ncat
57 CHARACTER ctype*4,phrase*100
60 REAL array(1000),param(3),prob(mxncat)
61 INTEGER iarray(1000),perm(500)
74 9000
FORMAT (
' Tests most generators of specific distributions.'/
75 +
' Generates 1000 deviates: reports mean and variance.'/
76 +
' Also reports theoretical mean and variance.'/
77 +
' If theoretical mean or var doesn''t exist prints -1.'/
78 +
' For permutations, generates one permutation of 1..n'/
80 +
' For uniform integers asks for upper bound, number of'/
81 +
' replicates per integer in 1..upper bound.'/
82 +
' Prints table of num times each integer generated.'/
83 +
' For multinomial asks for number of events to be'/
84 +
' classified, number of categories in which they'/
85 +
' are to be classified, and the probabilities that'/
86 +
' an event will be classified in the categories,'/
87 +
' for all but the last category. Prints table of'/
88 +
' number of events by category, true probability'/
89 +
' associated with each category, and observed'/
90 +
' proportion of events in each category.')
96 9010
FORMAT (
' Enter number corresponding to choice:'/
97 +
' (0) Exit this program'/
98 +
' (1) Generate Chi-Square deviates'/
99 +
' (2) Generate noncentral Chi-Square deviates'/
100 +
' (3) Generate F deviates'/
101 +
' (4) Generate noncentral F deviates'/
102 +
' (5) Generate random permutation'/
103 +
' (6) Generate uniform integers'/
104 +
' (7) Generate uniform reals'/
105 +
' (8) Generate beta deviates'/
106 +
' (9) Generate binomial outcomes'/
107 +
' (10) Generate Poisson outcomes'/
108 +
' (11) Generate exponential deviates'/
109 +
' (12) Generate gamma deviates'/
110 +
' (13) Generate multinomial outcomes'/
111 +
' (14) Generate normal deviates'/
112 +
' (15) Generate negative binomial outcomes'/)
115 IF (.NOT. (iwhich.LT.0.OR.iwhich.GT.mxwh))
GO TO 20
116 WRITE (*,*)
' Choices are 1..',mxwh,
' - try again.' 119 20
IF (iwhich.EQ.0) stop
' Normal termination rn tests' 120 WRITE (*,*)
' Enter phrase to initialize rn generator' 121 READ (*,
'(a)') phrase
122 CALL phrtsd(phrase,is1,is2)
125 IF ((1).NE. (iwhich))
GO TO 40
130 WRITE (*,*)
' Enter (real) df for the chi-square generation' 136 CALL trstat(ctype,param,avtr,vartr)
137 WRITE (*,9020) av,avtr,var,vartr,
xmin,
xmax 139 9020
FORMAT (
' Mean Generated: ',t30,g15.7,5x,
'True:',t60,
140 + g15.7/
' Variance Generated:',t30,g15.7,5x,
'True:',t60,
141 + g15.7/
' Minimum: ',t30,g15.7,5x,
'Maximum:',t60,g15.7)
145 40
IF ((2).NE. (iwhich))
GO TO 60
151 WRITE (*,*)
' Enter (real) df' 152 WRITE (*,*)
' (real) noncentrality parameter' 153 READ (*,*) param(1),param(2)
158 CALL trstat(ctype,param,avtr,vartr)
159 WRITE (*,9020) av,avtr,var,vartr,
xmin,
xmax 162 60
IF ((3).NE. (iwhich))
GO TO 80
168 WRITE (*,*)
' Enter (real) df of the numerator' 169 WRITE (*,*)
' (real) df of the denominator' 170 READ (*,*) param(1),param(2)
175 CALL trstat(ctype,param,avtr,vartr)
176 WRITE (*,9020) av,avtr,var,vartr,
xmin,
xmax 179 80
IF ((4).NE. (iwhich))
GO TO 100
185 WRITE (*,*)
' Enter (real) df of the numerator' 186 WRITE (*,*)
' (real) df of the denominator' 187 WRITE (*,*)
' (real) noncentrality parameter' 188 READ (*,*) param(1),param(2),param(3)
193 CALL trstat(ctype,param,avtr,vartr)
194 WRITE (*,9020) av,avtr,var,vartr,
xmin,
xmax 197 100
IF ((5).NE. (iwhich))
GO TO 140
202 110
WRITE (*,*)
' Enter size of permutation' 204 IF (.NOT. (nperm.LT.1.OR.nperm.GT.500))
GO TO 120
205 WRITE (*,*)
' Permutation size must be between 1 and 500 ',
209 120
WRITE (*,*)
' Random Permutation Generated - Size',nperm
214 WRITE (*,*)
' Perm Generated' 215 WRITE (*,
'(20I4)') (perm(i),i=1,nperm)
218 140
IF ((6).NE. (iwhich))
GO TO 170
223 WRITE (*,*)
' Enter maximum uniform integer' 225 WRITE (*,*)
' Enter number of replications per integer' 233 iarray(itmp) = iarray(itmp) + 1
235 WRITE (*,*)
' Counts of Integers Generated' 236 WRITE (*,
'(20I4)') (iarray(j),j=1,mxint)
239 170
IF ((7).NE. (iwhich))
GO TO 190
245 WRITE (*,*)
' Enter Low then High bound for uniforms' 246 READ (*,*) param(1),param(2)
251 CALL trstat(ctype,param,avtr,vartr)
252 WRITE (*,9020) av,avtr,var,vartr,
xmin,
xmax 255 190
IF ((8).NE. (iwhich))
GO TO 210
261 WRITE (*,*)
' Enter A, B for Beta deviate' 262 READ (*,*) param(1),param(2)
267 CALL trstat(ctype,param,avtr,vartr)
268 WRITE (*,9020) av,avtr,var,vartr,
xmin,
xmax 271 210
IF ((9).NE. (iwhich))
GO TO 240
277 WRITE (*,*)
' Enter number of trials, Prob event for ',
278 +
'binomial outcomes' 281 iarray(i) =
ignbin(ntry,pevt)
289 CALL trstat(ctype,param,avtr,vartr)
290 WRITE (*,9020) av,avtr,var,vartr,
xmin,
xmax 293 240
IF ((10).NE. (iwhich))
GO TO 270
299 WRITE (*,*)
' Enter mean for Poisson generation' 302 iarray(i) =
ignpoi(param(1))
308 CALL trstat(ctype,param,avtr,vartr)
309 WRITE (*,9020) av,avtr,var,vartr,
xmin,
xmax 312 270
IF ((11).NE. (iwhich))
GO TO 290
318 WRITE (*,*)
' Enter (real) AV for Exponential' 324 CALL trstat(ctype,param,avtr,vartr)
325 WRITE (*,9020) av,avtr,var,vartr,
xmin,
xmax 329 290
IF ((12).NE. (iwhich))
GO TO 310
335 WRITE (*,*)
' Enter (real) A, (real) R for Gamma deviate' 336 READ (*,*) param(1),param(2)
341 CALL trstat(ctype,param,avtr,vartr)
342 WRITE (*,9020) av,avtr,var,vartr,
xmin,
xmax 345 310
IF ((13).NE. (iwhich))
GO TO 360
350 WRITE (*,*)
' Enter (int) number of observations: ' 352 320
WRITE (*,*)
' Enter (int) num. of categories: <= ',mxncat
354 IF (ncat.GT.mxncat)
THEN 355 WRITE (*,*)
' number of categories must be <= ',mxncat
356 WRITE (*,*)
' Try again ... ' 359 WRITE (*,*)
' Enter (real) prob. vector of length ',ncat-1
360 READ (*,*) (prob(i),i=1,ncat-1)
361 CALL genmul(ntry,prob,ncat,iarray)
366 ntot = ntot + iarray(i)
367 array(i) = iarray(i)/rtry
371 ntot = ntot + iarray(i)
377 psum = psum + prob(i)
379 prob(ncat) = 1.0 - psum
381 WRITE (*,*)
' Total number of observations: ',ntot
382 WRITE (*,*)
' Total observations by category: ' 383 WRITE (*,
'(10I8)') (iarray(i),i=1,ncat)
384 WRITE (*,*)
' True probabilities by category: ' 385 WRITE (*,
'(8F10.7)') (prob(i),i=1,ncat)
386 WRITE (*,*)
' Observed proportions by category: ' 387 WRITE (*,
'(8F10.7)') (
array(i),i=1,ncat)
390 360
IF ((14).NE. (iwhich))
GO TO 380
396 WRITE (*,*)
' Enter (real) AV, (real) SD for Normal' 397 READ (*,*) param(1),param(2)
402 CALL trstat(ctype,param,avtr,vartr)
403 WRITE (*,9020) av,avtr,var,vartr,
xmin,
xmax 406 380
IF ((15).NE. (iwhich))
GO TO 410
412 WRITE (*,*)
' Enter required (int) Number of events then ' 413 WRITE (*,*)
' (real) Prob of an event for negative binomial' 416 iarray(i) =
ignnbn(ntry,pevt)
424 CALL trstat(ctype,param,avtr,vartr)
425 WRITE (*,9020) av,avtr,var,vartr,
xmin,
xmax 432 SUBROUTINE trstat(ctype,parin,av,var)
433 IMPLICIT INTEGER (i-n),
REAL (a-h,o-p,r-z),
LOGICAL (q)
522 IF ((
'chis').NE. (ctype))
GO TO 10
527 10
IF ((
'ncch').NE. (ctype))
GO TO 20
528 a = parin(1) + parin(2)
534 20
IF ((
'f').NE. (ctype))
GO TO 70
535 IF (.NOT. (parin(2).LE.2.0001))
GO TO 30
539 30 av = parin(2)/ (parin(2)-2.0)
540 40
IF (.NOT. (parin(2).LE.4.0001))
GO TO 50
544 50 var = (2.0*parin(2)**2* (parin(1)+parin(2)-2.0))/
545 + (parin(1)* (parin(2)-2.0)**2* (parin(2)-4.0))
548 70
IF ((
'ncf').NE. (ctype))
GO TO 120
549 IF (.NOT. (parin(2).LE.2.0001))
GO TO 80
553 80 av = (parin(2)* (parin(1)+parin(3)))/ ((parin(2)-2.0)*parin(1))
554 90
IF (.NOT. (parin(2).LE.4.0001))
GO TO 100
558 100 a = (parin(1)+parin(3))**2 + (parin(1)+2.0*parin(3))*
560 b = (parin(2)-2.0)**2* (parin(2)-4.0)
561 var = 2.0* (parin(2)/parin(1))**2* (a/b)
564 120
IF ((
'unif').NE. (ctype))
GO TO 130
565 range = parin(2) - parin(1)
566 av = parin(1) + range/2.0
570 130
IF ((
'beta').NE. (ctype))
GO TO 140
571 av = parin(1)/ (parin(1)+parin(2))
572 var = (av*parin(2))/ ((parin(1)+parin(2))*
573 + (parin(1)+parin(2)+1.0))
576 140
IF ((
'bin').NE. (ctype))
GO TO 150
577 av = parin(1)*parin(2)
578 var = av* (1.0-parin(2))
581 150
IF ((
'pois').NE. (ctype))
GO TO 160
586 160
IF ((
'expo').NE. (ctype))
GO TO 170
591 170
IF ((
'gamm').NE. (ctype))
GO TO 180
592 av = parin(2) / parin(1)
596 180
IF ((
'norm').NE. (ctype))
GO TO 190
601 190
IF ((
'nbin').NE. (ctype))
GO TO 200
602 av = parin(1) * (1.0 - parin(2)) / parin(2)
606 200
WRITE (*,*)
'Unimplemented type ',ctype
607 stop
'Unimplemented type in TRSTAT' real function genf(dfn, dfd)
real function gennor(av, sd)
octave_int< T > xmax(const octave_int< T > &x, const octave_int< T > &y)
real function genunf(low, high)
then the function must return scalars which will be concatenated into the return array(s). If code
integer *4 function ignbin(n, pp)
subroutine phrtsd(phrase, seed1, seed2)
real function genbet(aa, bb)
integer *4 function ignnbn(n, p)
subroutine genprm(iarray, larray)
integer *4 function ignuin(low, high)
subroutine genmul(n, p, ncat, ix)
real function gengam(a, r)
real function gennch(df, xnonc)
subroutine setall(iseed1, iseed2)
octave_int< T > xmin(const octave_int< T > &x, const octave_int< T > &y)
real function gennf(dfn, dfd, xnonc)
integer *4 function ignpoi(mu)
subroutine stat(x, n, av, var, xmin, xmax)
ColumnVector real(const ComplexColumnVector &a)
subroutine trstat(ctype, parin, av, var)