GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
tstgmn.for
Go to the documentation of this file.
1 C JJV changed name to ONECOV to avoid confusion with array COVAR
2 C JJV this was also changed in the body of the function
3 C REAL FUNCTION covar(x,y,n)
4  REAL FUNCTION onecov(x,y,n)
5 C .. Scalar Arguments ..
6  INTEGER n
7 C ..
8 C .. Array Arguments ..
9  REAL x(n),y(n)
10 C ..
11 C .. Local Scalars ..
12  REAL avx,avy,varx,vary,xmax,xmin
13  INTEGER i
14 C ..
15 C .. External Subroutines ..
16  EXTERNAL stat
17 C ..
18 C .. Intrinsic Functions ..
19  INTRINSIC real
20 C ..
21 C .. Executable Statements ..
22  CALL stat(x,n,avx,varx,xmin,xmax)
23  CALL stat(y,n,avy,vary,xmin,xmax)
24 C covar = 0.0
25  onecov = 0.0
26  DO 10,i = 1,n
27 C covar = covar + (x(i)-avx)* (y(i)-avy)
28  onecov = onecov + (x(i)-avx)* (y(i)-avy)
29  10 CONTINUE
30 C covar = covar/real(n-1)
31  onecov = onecov/real(n-1)
32  RETURN
33 
34  END
35 
36 C JJV Added argument LDXCOV (leading dimension of XCOVAR) to be
37 C JJV consistent with the program TSTGMN, see comments below.
38 C JJV This change necessitated changes in the declarations.
39 C SUBROUTINE prcomp(p,mean,xcovar,answer)
40  SUBROUTINE prcomp(p,mean,xcovar,ldxcov,answer)
41 
42 C INTEGER p,maxp
43  INTEGER p,maxp,ldxcov
44  parameter(maxp=10)
45 C REAL mean(p),xcovar(p,p),rcovar(maxp,maxp)
46  REAL mean(p),xcovar(ldxcov,p),rcovar(maxp,maxp)
47  REAL answer(1000,maxp)
48 C JJV added ONECOV because of name change to function COVAR
49 C REAL rmean(maxp),rvar(maxp)
50  REAL rmean(maxp),rvar(maxp),onecov
51  INTEGER maxobs
52  parameter(maxobs=1000)
53 
54  DO 10,i = 1,p
55  CALL stat(answer(1,i),maxobs,rmean(i),rvar(i),dum1,dum2)
56  WRITE (*,*) ' Variable Number',i
57  WRITE (*,*) ' Mean ',mean(i),' Generated ',rmean(i)
58  WRITE (*,*) ' Variance ',xcovar(i,i),' Generated',rvar(i)
59  10 CONTINUE
60  WRITE (*,*) ' Covariances'
61  DO 30,i = 1,p
62  DO 20,j = 1,i - 1
63  WRITE (*,*) ' I = ',i,' J = ',j
64 C JJV changed COVAR to match new name
65 C rcovar(i,j) = covar(answer(1,i),answer(1,j),maxobs)
66  rcovar(i,j) = onecov(answer(1,i),answer(1,j),maxobs)
67  WRITE (*,*) ' Covariance ',xcovar(i,j),' Generated ',
68  + rcovar(i,j)
69  20 CONTINUE
70  30 CONTINUE
71  RETURN
72 
73  END
74 
75 C JJV added LDCOV (leading dimension of COVAR) to be
76 C JJV consistent with the program TSTGMN, see comments below.
77 C JJV This change necessitated changes in the declarations.
78 C SUBROUTINE setcov(p,var,corr,covar)
79  SUBROUTINE setcov(p,var,corr,covar,ldcov)
80 C Set covariance matrix from variance and common correlation
81 C .. Scalar Arguments ..
82  REAL corr
83 C INTEGER p
84  INTEGER p,ldcov
85 C ..
86 C .. Array Arguments ..
87 C REAL covar(p,p),var(p)
88  REAL covar(ldcov,p),var(p)
89 C ..
90 C .. Local Scalars ..
91  INTEGER i,j
92 C ..
93 C .. Intrinsic Functions ..
94  INTRINSIC sqrt
95 C ..
96 C .. Executable Statements ..
97  DO 40,i = 1,p
98  DO 30,j = 1,p
99  IF (.NOT. (i.EQ.j)) GO TO 10
100  covar(i,j) = var(i)
101  GO TO 20
102 
103  10 covar(i,j) = corr*sqrt(var(i)*var(j))
104  20 CONTINUE
105  30 CONTINUE
106  40 CONTINUE
107  RETURN
108 
109  END
110 
111  SUBROUTINE stat(x,n,av,var,xmin,xmax)
112 C .. Scalar Arguments ..
113  REAL av,var,xmax,xmin
114  INTEGER n
115 C ..
116 C .. Array Arguments ..
117  REAL x(n)
118 C ..
119 C .. Local Scalars ..
120  REAL sum
121  INTEGER i
122 C ..
123 C .. Intrinsic Functions ..
124  INTRINSIC real
125 C ..
126 C .. Executable Statements ..
127  xmin = x(1)
128  xmax = x(1)
129  sum = 0.0
130  DO 10,i = 1,n
131  sum = sum + x(i)
132  IF (x(i).LT.xmin) xmin = x(i)
133  IF (x(i).GT.xmax) xmax = x(i)
134  10 CONTINUE
135  av = sum/real(n)
136  sum = 0.0
137  DO 20,i = 1,n
138  sum = sum + (x(i)-av)**2
139  20 CONTINUE
140  var = sum/real(n-1)
141  RETURN
142 
143  END
144 
145  PROGRAM tstgmn
146 C Test Generation of Multivariate Normal Data
147 C JJV SETGMN was: SUBROUTINE setgmn(meanv,covm,p,parm)
148 C JJV is: SUBROUTINE setgmn(meanv,covm,ldcovm,p,parm)
149 C JJV So the covariance matrices have been changed to 2-dim'l
150 C JJV matrices, and the additional argument has been added to
151 C JJV the subroutine call. Additional changes have been made
152 C JJV to reflect this. (in declarations, the matrix copy routine,
153 C JJV and in subroutine calls.)
154 C .. Parameters ..
155  INTEGER maxp
156  parameter(maxp=10)
157  INTEGER maxobs
158  parameter(maxobs=1000)
159 C JJV this parameter is no longer needed
160 C INTEGER p2
161 C PARAMETER (p2=maxp*maxp)
162 C ..
163 C .. Local Scalars ..
164  REAL corr
165  INTEGER i,iobs,is1,is2,j,p
166  CHARACTER phrase*100
167 C ..
168 C .. Local Arrays ..
169 C REAL answer(1000,maxp),ccovar(p2),covar(p2),mean(maxp),param(500),
170 C + temp(maxp),var(maxp),work(maxp)
171  REAL answer(1000,maxp),ccovar(maxp,maxp),covar(maxp,maxp),
172  + mean(maxp),param(500),temp(maxp),var(maxp),work(maxp)
173 C ..
174 C .. External Subroutines ..
175  EXTERNAL genmn,phrtsd,prcomp,setall,setcov,setgmn
176 C ..
177 C .. Executable Statements ..
178  WRITE (*,9000)
179 
180  9000 FORMAT (
181  + ' Tests Multivariate Normal Generator for Up to 10 Variables'
182  + /
183  + ' User inputs means, variances, one correlation that is applied'
184  + /' to all pairs of variables'/
185  + ' 1000 multivariate normal deviates are generated'/
186  + ' Means, variances and covariances are calculated for these.'
187  + )
188 
189  10 WRITE (*,*) 'Enter number of variables for normal generator'
190  READ (*,*) p
191  WRITE (*,*) 'Enter mean vector of length ',p
192  READ (*,*) (mean(i),i=1,p)
193  WRITE (*,*) 'Enter variance vector of length ',p
194  READ (*,*) (var(i),i=1,p)
195  WRITE (*,*) 'Enter correlation of all variables'
196  READ (*,*) corr
197 C CALL setcov(p,var,corr,covar)
198  CALL setcov(p,var,corr,covar,maxp)
199  WRITE (*,*) ' Enter phrase to initialize rn generator'
200  READ (*,'(a)') phrase
201  CALL phrtsd(phrase,is1,is2)
202  CALL setall(is1,is2)
203 C DO 20,i = 1,p2
204 C ccovar(i) = covar(i)
205 C 20 CONTINUE
206  DO 25,i = 1,maxp
207  DO 20,j = 1,maxp
208  ccovar(i,j) = covar(i,j)
209  20 CONTINUE
210  25 CONTINUE
211 C
212 C Generate Variables
213 C
214 C CALL setgmn(mean,ccovar,p,param)
215  CALL setgmn(mean,ccovar,maxp,p,param)
216  DO 40,iobs = 1,maxobs
217  CALL genmn(param,work,temp)
218  DO 30,j = 1,p
219  answer(iobs,j) = work(j)
220  30 CONTINUE
221  40 CONTINUE
222 C CALL prcomp(p,mean,covar,answer)
223  CALL prcomp(p,mean,covar,maxp,answer)
224 C
225 C Print Comparison of Generated and Reconstructed Values
226 C
227  GO TO 10
228 
229  END
xmax
octave_int< T > xmax(const octave_int< T > &x, const octave_int< T > &y)
Definition: oct-inttypes.h:1224
setcov
subroutine setcov(p, var, corr, covar, ldcov)
Definition: tstgmn.for:80
phrtsd
subroutine phrtsd(phrase, seed1, seed2)
Definition: phrtsd.f:2
onecov
real function onecov(x, y, n)
Definition: tstgmn.for:5
tstgmn
program tstgmn
Definition: tstgmn.for:145
prcomp
subroutine prcomp(p, mean, xcovar, ldxcov, answer)
Definition: tstgmn.for:41
setgmn
subroutine setgmn(meanv, covm, ldcovm, p, parm)
Definition: setgmn.f:2
setall
subroutine setall(iseed1, iseed2)
Definition: setall.f:2
xmin
octave_int< T > xmin(const octave_int< T > &x, const octave_int< T > &y)
Definition: oct-inttypes.h:1233
stat
subroutine stat(x, n, av, var, xmin, xmax)
Definition: tstgmn.for:112
genmn
subroutine genmn(parm, x, work)
Definition: genmn.f:2
real
ColumnVector real(const ComplexColumnVector &a)
Definition: dColVector.cc:135