|
2330
|
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 .. |
|
3188
|
51 INTEGER mxwh,mxncat |
|
|
52 PARAMETER (mxwh=15,mxncat=100) |
|
2330
|
53 C .. |
|
|
54 C .. Local Scalars .. |
|
3188
|
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 |
|
2330
|
57 CHARACTER type*4,phrase*100 |
|
|
58 C .. |
|
|
59 C .. Local Arrays .. |
|
3188
|
60 REAL array(1000),param(3),prob(mxncat) |
|
2330
|
61 INTEGER iarray(1000),perm(500) |
|
|
62 C .. |
|
|
63 C .. External Functions .. |
|
3188
|
64 REAL genbet,genchi,genf,gennch,gennf,genunf,genexp,gengam,gennor |
|
|
65 INTEGER ignuin,ignnbn |
|
|
66 EXTERNAL genbet,genchi,genf,gennch,gennf,genunf,ignuin |
|
2330
|
67 C .. |
|
|
68 C .. External Subroutines .. |
|
3188
|
69 EXTERNAL genprm,phrtsd,setall,stat,trstat,genmul |
|
2330
|
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.'/ |
|
3188
|
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.') |
|
2330
|
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'/ |
|
3188
|
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'/) |
|
2330
|
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, |
|
3188
|
141 + G15.7/' Minimum: ',T30,G15.7,5X,'Maximum:',T60,G15.7) |
|
2330
|
142 |
|
3188
|
143 GO TO 420 |
|
2330
|
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 |
|
3188
|
160 GO TO 420 |
|
2330
|
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 |
|
3188
|
177 GO TO 420 |
|
2330
|
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 |
|
3188
|
195 GO TO 420 |
|
2330
|
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) |
|
3188
|
216 GO TO 420 |
|
2330
|
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) |
|
3188
|
237 GO TO 420 |
|
2330
|
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 |
|
3188
|
253 GO TO 420 |
|
2330
|
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 |
|
3188
|
269 GO TO 420 |
|
2330
|
270 |
|
|
271 210 IF ((9).NE. (iwhich)) GO TO 240 |
|
3188
|
272 |
|
2330
|
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 |
|
3188
|
291 GO TO 420 |
|
2330
|
292 |
|
|
293 240 IF ((10).NE. (iwhich)) GO TO 270 |
|
3188
|
294 |
|
2330
|
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 |
|
3188
|
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 |
|
2330
|
346 |
|
3188
|
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 |
|
2330
|
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 |
|
3188
|
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 |
|
2330
|
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 |
|
3188
|
480 C binomial |
|
2330
|
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 |
|
3188
|
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(*) |
|
2330
|
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) |
|
3188
|
525 GO TO 210 |
|
2330
|
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) |
|
3188
|
532 GO TO 210 |
|
2330
|
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)) |
|
3188
|
546 60 GO TO 210 |
|
2330
|
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) |
|
3188
|
562 110 GO TO 210 |
|
2330
|
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 |
|
3188
|
568 GO TO 210 |
|
2330
|
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)) |
|
3188
|
574 GO TO 210 |
|
2330
|
575 |
|
|
576 140 IF (('bin').NE. (type)) GO TO 150 |
|
|
577 av = parin(1)*parin(2) |
|
|
578 var = av* (1.0-parin(2)) |
|
3188
|
579 GO TO 210 |
|
2330
|
580 |
|
|
581 150 IF (('pois').NE. (type)) GO TO 160 |
|
|
582 av = parin(1) |
|
|
583 var = parin(1) |
|
3188
|
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 |
|
2330
|
595 |
|
3188
|
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 |
|
2330
|
607 STOP 'Unimplemented type in TRSTAT' |
|
|
608 |
|
3188
|
609 210 RETURN |
|
2330
|
610 |
|
|
611 END |
