GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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 type*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  type = '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(type,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  type = '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(type,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  type = '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(type,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  type = '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(type,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  type = '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(type,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  type = '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(type,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  type = '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(type,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  type = '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(type,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  type = '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(type,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  type = '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(type,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  type = '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(type,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  type = '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(type,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(type,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 TYPE --> 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 type* (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. (type)) go to 10
523  av = parin(1)
524  var = 2.0*parin(1)
525  go to 210
526 
527  10 IF (('ncch').NE. (type)) 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. (type)) 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. (type)) 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. (type)) 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. (type)) 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. (type)) 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. (type)) go to 160
582  av = parin(1)
583  var = parin(1)
584  go to 210
585 
586  160 IF (('expo').NE. (type)) go to 170
587  av = parin(1)
588  var = parin(1)**2
589  go to 210
590 
591  170 IF (('gamm').NE. (type)) go to 180
592  av = parin(2) / parin(1)
593  var = av / parin(1)
594  go to 210
595 
596  180 IF (('norm').NE. (type)) go to 190
597  av = parin(1)
598  var = parin(2)**2
599  go to 210
600 
601  190 IF (('nbin').NE. (type)) 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 ',type
607  stop 'Unimplemented type in TRSTAT'
608 
609  210 RETURN
610 
611  END
Generated on Mon Dec 30 2013 03:04:47 for GNU Octave by   doxygen 1.8.1.2