GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
tstmid.for
Go to the documentation of this file.
1  SUBROUTINE stat(x,n,av,var,xmin,xmax)
2 C**********************************************************************
3 C
4 C SUBROUTINE STAT( X, N, AV, VAR)
5 C
6 C compute STATistics
7 C
8 C
9 C Function
10 C
11 C
12 C Computes AVerage and VARiance of array X(N).
13 C
14 C**********************************************************************
15 C .. Scalar Arguments ..
16  REAL av,var,xmax,xmin
17  INTEGER n
18 C ..
19 C .. Array Arguments ..
20  REAL x(n)
21 C ..
22 C .. Local Scalars ..
23  REAL sum
24  INTEGER i
25 C ..
26 C .. Intrinsic Functions ..
27  INTRINSIC real
28 C ..
29 C .. Executable Statements ..
30  xmin = x(1)
31  xmax = x(1)
32  sum = 0.0
33  DO 10,i = 1,n
34  sum = sum + x(i)
35  IF (x(i).LT.xmin) xmin = x(i)
36  IF (x(i).GT.xmax) xmax = x(i)
37  10 CONTINUE
38  av = sum/real(n)
39  sum = 0.0
40  DO 20,i = 1,n
41  sum = sum + (x(i)-av)**2
42  20 CONTINUE
43  var = sum/real(n-1)
44  RETURN
45 
46  END
47  PROGRAM tstall
48  IMPLICIT LOGICAL (q)
49 C Interactive test for PHRTSD
50 C .. Parameters ..
51  INTEGER mxwh,mxncat
52  parameter(mxwh=15,mxncat=100)
53 C ..
54 C .. Local Scalars ..
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
58 C ..
59 C .. Local Arrays ..
60  REAL array(1000),param(3),prob(mxncat)
61  INTEGER iarray(1000),perm(500)
62 C ..
63 C .. External Functions ..
64  REAL genbet,genchi,genf,gennch,gennf,genunf,genexp,gengam,gennor
65  INTEGER ignuin,ignnbn
66  EXTERNAL genbet,genchi,genf,gennch,gennf,genunf,ignuin
67 C ..
68 C .. External Subroutines ..
69  EXTERNAL genprm,phrtsd,setall,stat,trstat,genmul
70 C ..
71 C .. Executable Statements ..
72  WRITE (*,9000)
73 
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'/
79  + ' and prints it.'/
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.')
91 C
92 C Menu for choosing tests
93 C
94  10 WRITE (*,9010)
95 
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'/)
113 
114  READ (*,*) iwhich
115  IF (.NOT. (iwhich.LT.0.OR.iwhich.GT.mxwh)) GO TO 20
116  WRITE (*,*) ' Choices are 1..',mxwh,' - try again.'
117  GO TO 10
118 
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)
123  CALL setall(is1,is2)
124 
125  IF ((1).NE. (iwhich)) GO TO 40
126 C
127 C Chi-square deviates
128 C
129  ctype = 'chis'
130  WRITE (*,*) ' Enter (real) df for the chi-square generation'
131  READ (*,*) param(1)
132  DO 30,i = 1,1000
133  array(i) = genchi(param(1))
134  30 CONTINUE
135  CALL stat(array,1000,av,var,xmin,xmax)
136  CALL trstat(ctype,param,avtr,vartr)
137  WRITE (*,9020) av,avtr,var,vartr,xmin,xmax
138 
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)
142 
143  GO TO 420
144 
145  40 IF ((2).NE. (iwhich)) GO TO 60
146 
147 C
148 C Noncentral Chi-square deviates
149 C
150  ctype = 'ncch'
151  WRITE (*,*) ' Enter (real) df'
152  WRITE (*,*) ' (real) noncentrality parameter'
153  READ (*,*) param(1),param(2)
154  DO 50,i = 1,1000
155  array(i) = gennch(param(1),param(2))
156  50 CONTINUE
157  CALL stat(array,1000,av,var,xmin,xmax)
158  CALL trstat(ctype,param,avtr,vartr)
159  WRITE (*,9020) av,avtr,var,vartr,xmin,xmax
160  GO TO 420
161 
162  60 IF ((3).NE. (iwhich)) GO TO 80
163 
164 C
165 C F deviates
166 C
167  ctype = 'f'
168  WRITE (*,*) ' Enter (real) df of the numerator'
169  WRITE (*,*) ' (real) df of the denominator'
170  READ (*,*) param(1),param(2)
171  DO 70,i = 1,1000
172  array(i) = genf(param(1),param(2))
173  70 CONTINUE
174  CALL stat(array,1000,av,var,xmin,xmax)
175  CALL trstat(ctype,param,avtr,vartr)
176  WRITE (*,9020) av,avtr,var,vartr,xmin,xmax
177  GO TO 420
178 
179  80 IF ((4).NE. (iwhich)) GO TO 100
180 
181 C
182 C Noncentral F deviates
183 C
184  ctype = 'ncf'
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)
189  DO 90,i = 1,1000
190  array(i) = gennf(param(1),param(2),param(3))
191  90 CONTINUE
192  CALL stat(array,1000,av,var,xmin,xmax)
193  CALL trstat(ctype,param,avtr,vartr)
194  WRITE (*,9020) av,avtr,var,vartr,xmin,xmax
195  GO TO 420
196 
197  100 IF ((5).NE. (iwhich)) GO TO 140
198 
199 C
200 C Random permutation
201 C
202  110 WRITE (*,*) ' Enter size of permutation'
203  READ (*,*) nperm
204  IF (.NOT. (nperm.LT.1.OR.nperm.GT.500)) GO TO 120
205  WRITE (*,*) ' Permutation size must be between 1 and 500 ',
206  + '- try again!'
207  GO TO 110
208 
209  120 WRITE (*,*) ' Random Permutation Generated - Size',nperm
210  DO 130,i = 1,500
211  perm(i) = i
212  130 CONTINUE
213  CALL genprm(perm,nperm)
214  WRITE (*,*) ' Perm Generated'
215  WRITE (*,'(20I4)') (perm(i),i=1,nperm)
216  GO TO 420
217 
218  140 IF ((6).NE. (iwhich)) GO TO 170
219 
220 C
221 C Uniform integer
222 C
223  WRITE (*,*) ' Enter maximum uniform integer'
224  READ (*,*) mxint
225  WRITE (*,*) ' Enter number of replications per integer'
226  READ (*,*) nrep
227  DO 150,i = 1,1000
228  iarray(i) = 0
229  150 CONTINUE
230  ntot = mxint*nrep
231  DO 160,i = 1,ntot
232  itmp = ignuin(1,mxint)
233  iarray(itmp) = iarray(itmp) + 1
234  160 CONTINUE
235  WRITE (*,*) ' Counts of Integers Generated'
236  WRITE (*,'(20I4)') (iarray(j),j=1,mxint)
237  GO TO 420
238 
239  170 IF ((7).NE. (iwhich)) GO TO 190
240 
241 C
242 C Uniform real
243 C
244  ctype = 'unif'
245  WRITE (*,*) ' Enter Low then High bound for uniforms'
246  READ (*,*) param(1),param(2)
247  DO 180,i = 1,1000
248  array(i) = genunf(param(1),param(2))
249  180 CONTINUE
250  CALL stat(array,1000,av,var,xmin,xmax)
251  CALL trstat(ctype,param,avtr,vartr)
252  WRITE (*,9020) av,avtr,var,vartr,xmin,xmax
253  GO TO 420
254 
255  190 IF ((8).NE. (iwhich)) GO TO 210
256 
257 C
258 C Beta deviate
259 C
260  ctype = 'beta'
261  WRITE (*,*) ' Enter A, B for Beta deviate'
262  READ (*,*) param(1),param(2)
263  DO 200,i = 1,1000
264  array(i) = genbet(param(1),param(2))
265  200 CONTINUE
266  CALL stat(array,1000,av,var,xmin,xmax)
267  CALL trstat(ctype,param,avtr,vartr)
268  WRITE (*,9020) av,avtr,var,vartr,xmin,xmax
269  GO TO 420
270 
271  210 IF ((9).NE. (iwhich)) GO TO 240
272 
273 C
274 C Binomial outcomes
275 C
276  ctype = 'bin'
277  WRITE (*,*) ' Enter number of trials, Prob event for ',
278  + 'binomial outcomes'
279  READ (*,*) ntry,pevt
280  DO 220,i = 1,1000
281  iarray(i) = ignbin(ntry,pevt)
282  220 CONTINUE
283  DO 230,i = 1,1000
284  array(i) = iarray(i)
285  230 CONTINUE
286  CALL stat(array,1000,av,var,xmin,xmax)
287  param(1) = ntry
288  param(2) = pevt
289  CALL trstat(ctype,param,avtr,vartr)
290  WRITE (*,9020) av,avtr,var,vartr,xmin,xmax
291  GO TO 420
292 
293  240 IF ((10).NE. (iwhich)) GO TO 270
294 
295 C
296 C Poisson outcomes
297 C
298  ctype = 'pois'
299  WRITE (*,*) ' Enter mean for Poisson generation'
300  READ (*,*) param(1)
301  DO 250,i = 1,1000
302  iarray(i) = ignpoi(param(1))
303  250 CONTINUE
304  DO 260,i = 1,1000
305  array(i) = iarray(i)
306  260 CONTINUE
307  CALL stat(array,1000,av,var,xmin,xmax)
308  CALL trstat(ctype,param,avtr,vartr)
309  WRITE (*,9020) av,avtr,var,vartr,xmin,xmax
310  GO TO 420
311 
312  270 IF ((11).NE. (iwhich)) GO TO 290
313 
314 C
315 C Exponential deviates
316 C
317  ctype = 'expo'
318  WRITE (*,*) ' Enter (real) AV for Exponential'
319  READ (*,*) param(1)
320  DO 280,i = 1,1000
321  array(i) = genexp(param(1))
322  280 CONTINUE
323  CALL stat(array,1000,av,var,xmin,xmax)
324  CALL trstat(ctype,param,avtr,vartr)
325  WRITE (*,9020) av,avtr,var,vartr,xmin,xmax
326 
327  GO TO 420
328 
329  290 IF ((12).NE. (iwhich)) GO TO 310
330 
331 C
332 C Gamma deviates
333 C
334  ctype = 'gamm'
335  WRITE (*,*) ' Enter (real) A, (real) R for Gamma deviate'
336  READ (*,*) param(1),param(2)
337  DO 300,i = 1,1000
338  array(i) = gengam(param(1),param(2))
339  300 CONTINUE
340  CALL stat(array,1000,av,var,xmin,xmax)
341  CALL trstat(ctype,param,avtr,vartr)
342  WRITE (*,9020) av,avtr,var,vartr,xmin,xmax
343  GO TO 420
344 
345  310 IF ((13).NE. (iwhich)) GO TO 360
346 
347 C
348 C Multinomial outcomes
349 C
350  WRITE (*,*) ' Enter (int) number of observations: '
351  READ (*,*) ntry
352  320 WRITE (*,*) ' Enter (int) num. of categories: <= ',mxncat
353  READ (*,*) ncat
354  IF (ncat.GT.mxncat) THEN
355  WRITE (*,*) ' number of categories must be <= ',mxncat
356  WRITE (*,*) ' Try again ... '
357  GO TO 320
358  END IF
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)
362  ntot = 0
363  IF (ntry.GT.0) THEN
364  rtry = real(ntry)
365  DO 330, i = 1,ncat
366  ntot = ntot + iarray(i)
367  array(i) = iarray(i)/rtry
368  330 CONTINUE
369  ELSE
370  DO 340, i = 1,ncat
371  ntot = ntot + iarray(i)
372  array(i) = 0.0
373  340 CONTINUE
374  ENDIF
375  psum = 0.0
376  DO 350, i = 1,ncat-1
377  psum = psum + prob(i)
378  350 CONTINUE
379  prob(ncat) = 1.0 - psum
380 
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)
388  GO TO 420
389 
390  360 IF ((14).NE. (iwhich)) GO TO 380
391 
392 C
393 C Normal deviates
394 C
395  ctype = 'norm'
396  WRITE (*,*) ' Enter (real) AV, (real) SD for Normal'
397  READ (*,*) param(1),param(2)
398  DO 370,i = 1,1000
399  array(i) = gennor(param(1),param(2))
400  370 CONTINUE
401  CALL stat(array,1000,av,var,xmin,xmax)
402  CALL trstat(ctype,param,avtr,vartr)
403  WRITE (*,9020) av,avtr,var,vartr,xmin,xmax
404  GO TO 420
405 
406  380 IF ((15).NE. (iwhich)) GO TO 410
407 
408 C
409 C Negative Binomial outcomes
410 C
411  ctype = 'nbin'
412  WRITE (*,*) ' Enter required (int) Number of events then '
413  WRITE (*,*) ' (real) Prob of an event for negative binomial'
414  READ (*,*) ntry,pevt
415  DO 390,i = 1,1000
416  iarray(i) = ignnbn(ntry,pevt)
417  390 CONTINUE
418  DO 400,i = 1,1000
419  array(i) = iarray(i)
420  400 CONTINUE
421  CALL stat(array,1000,av,var,xmin,xmax)
422  param(1) = ntry
423  param(2) = pevt
424  CALL trstat(ctype,param,avtr,vartr)
425  WRITE (*,9020) av,avtr,var,vartr,xmin,xmax
426  GO TO 420
427 
428  410 CONTINUE
429  420 GO TO 10
430 
431  END
432  SUBROUTINE trstat(ctype,parin,av,var)
433  IMPLICIT INTEGER (i-n),real(a-h,o-p,r-z),logical(q)
434 C**********************************************************************
435 C
436 C SUBROUTINE TRSTAT( TYPE, PARIN, AV, VAR )
437 C TRue STATistics
438 C
439 C Returns mean and variance for a number of statistical distribution
440 C as a function of their parameters.
441 C
442 C
443 C Arguments
444 C
445 C
446 C CTYPE --> Character string indicating type of distribution
447 C 'chis' chisquare
448 C 'ncch' noncentral chisquare
449 C 'f' F (variance ratio)
450 C 'ncf' noncentral f
451 C 'unif' uniform
452 C 'beta' beta distribution
453 C 'bin' binomial
454 C 'pois' poisson
455 C 'expo' exponential
456 C 'gamm' gamma
457 C 'norm' normal
458 C 'nbin' negative binomial
459 C CHARACTER*(4) TYPE
460 C
461 C PARIN --> Array containing parameters of distribution
462 C chisquare
463 C PARIN(1) is df
464 C noncentral chisquare
465 C PARIN(1) is df
466 C PARIN(2) is noncentrality parameter
467 C F (variance ratio)
468 C PARIN(1) is df numerator
469 C PARIN(2) is df denominator
470 C noncentral F
471 C PARIN(1) is df numerator
472 C PARIN(2) is df denominator
473 C PARIN(3) is noncentrality parameter
474 C uniform
475 C PARIN(1) is LOW bound
476 C PARIN(2) is HIGH bound
477 C beta
478 C PARIN(1) is A
479 C PARIN(2) is B
480 C binomial
481 C PARIN(1) is Number of trials
482 C PARIN(2) is Prob Event at Each Trial
483 C poisson
484 C PARIN(1) is Mean
485 C exponential
486 C PARIN(1) is Mean
487 C gamma
488 C PARIN(1) is A
489 C PARIN(2) is R
490 C normal
491 C PARIN(1) is Mean
492 C PARIN(2) is Standard Deviation
493 C negative binomial
494 C PARIN(1) is required Number of events
495 C PARIN(2) is Probability of event
496 C REAL PARIN(*)
497 C
498 C AV <-- Mean of specified distribution with specified parameters
499 C REAL AV
500 C
501 C VAR <-- Variance of specified distribution with specified paramete
502 C REAL VAR
503 C
504 C
505 C Note
506 C
507 C
508 C AV and Var will be returned -1 if mean or variance is infinite
509 C
510 C**********************************************************************
511 C .. Scalar Arguments ..
512  REAL av,var
513  CHARACTER ctype* (4)
514 C ..
515 C .. Array Arguments ..
516  REAL parin(*)
517 C ..
518 C .. Local Scalars ..
519  REAL a,b,range
520 C ..
521 C .. Executable Statements ..
522  IF (('chis').NE. (ctype)) GO TO 10
523  av = parin(1)
524  var = 2.0*parin(1)
525  GO TO 210
526 
527  10 IF (('ncch').NE. (ctype)) GO TO 20
528  a = parin(1) + parin(2)
529  b = parin(2)/a
530  av = a
531  var = 2.0*a* (1.0+b)
532  GO TO 210
533 
534  20 IF (('f').NE. (ctype)) GO TO 70
535  IF (.NOT. (parin(2).LE.2.0001)) GO TO 30
536  av = -1.0
537  GO TO 40
538 
539  30 av = parin(2)/ (parin(2)-2.0)
540  40 IF (.NOT. (parin(2).LE.4.0001)) GO TO 50
541  var = -1.0
542  GO TO 60
543 
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))
546  60 GO TO 210
547 
548  70 IF (('ncf').NE. (ctype)) GO TO 120
549  IF (.NOT. (parin(2).LE.2.0001)) GO TO 80
550  av = -1.0
551  GO TO 90
552 
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
555  var = -1.0
556  GO TO 110
557 
558  100 a = (parin(1)+parin(3))**2 + (parin(1)+2.0*parin(3))*
559  + (parin(2)-2.0)
560  b = (parin(2)-2.0)**2* (parin(2)-4.0)
561  var = 2.0* (parin(2)/parin(1))**2* (a/b)
562  110 GO TO 210
563 
564  120 IF (('unif').NE. (ctype)) GO TO 130
565  range = parin(2) - parin(1)
566  av = parin(1) + range/2.0
567  var = range**2/12.0
568  GO TO 210
569 
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))
574  GO TO 210
575 
576  140 IF (('bin').NE. (ctype)) GO TO 150
577  av = parin(1)*parin(2)
578  var = av* (1.0-parin(2))
579  GO TO 210
580 
581  150 IF (('pois').NE. (ctype)) GO TO 160
582  av = parin(1)
583  var = parin(1)
584  GO TO 210
585 
586  160 IF (('expo').NE. (ctype)) GO TO 170
587  av = parin(1)
588  var = parin(1)**2
589  GO TO 210
590 
591  170 IF (('gamm').NE. (ctype)) GO TO 180
592  av = parin(2) / parin(1)
593  var = av / parin(1)
594  GO TO 210
595 
596  180 IF (('norm').NE. (ctype)) GO TO 190
597  av = parin(1)
598  var = parin(2)**2
599  GO TO 210
600 
601  190 IF (('nbin').NE. (ctype)) GO TO 200
602  av = parin(1) * (1.0 - parin(2)) / parin(2)
603  var = av / parin(2)
604  GO TO 210
605 
606  200 WRITE (*,*) 'Unimplemented type ',ctype
607  stop 'Unimplemented type in TRSTAT'
608 
609  210 RETURN
610 
611  END
real
ColumnVector real(const ComplexColumnVector &a)
Definition: dColVector.cc:137
genbet
real function genbet(aa, bb)
Definition: genbet.f:2
genchi
real function genchi(df)
Definition: genchi.f:2
genexp
real function genexp(av)
Definition: genexp.f:2
genf
real function genf(dfn, dfd)
Definition: genf.f:2
gengam
real function gengam(a, r)
Definition: gengam.f:2
genmul
subroutine genmul(n, p, ncat, ix)
Definition: genmul.f:2
gennch
real function gennch(df, xnonc)
Definition: gennch.f:2
gennf
real function gennf(dfn, dfd, xnonc)
Definition: gennf.f:2
gennor
real function gennor(av, sd)
Definition: gennor.f:2
genprm
subroutine genprm(iarray, larray)
Definition: genprm.f:2
genunf
real function genunf(low, high)
Definition: genunf.f:2
ignbin
integer *4 function ignbin(n, pp)
Definition: ignbin.f:2
ignnbn
integer *4 function ignnbn(n, p)
Definition: ignnbn.f:2
ignpoi
integer *4 function ignpoi(mu)
Definition: ignpoi.f:2
ignuin
integer *4 function ignuin(low, high)
Definition: ignuin.f:2
xmin
octave_int< T > xmin(const octave_int< T > &x, const octave_int< T > &y)
Definition: oct-inttypes.h:1318
xmax
octave_int< T > xmax(const octave_int< T > &x, const octave_int< T > &y)
Definition: oct-inttypes.h:1308
phrtsd
subroutine phrtsd(phrase, seed1, seed2)
Definition: phrtsd.f:2
setall
subroutine setall(iseed1, iseed2)
Definition: setall.f:2
stat
subroutine stat(x, n, av, var, xmin, xmax)
Definition: tstgmn.for:112
tstall
program tstall
Definition: tstmid.for:47
trstat
subroutine trstat(ctype, parin, av, var)
Definition: tstmid.for:433