GNU Octave 10.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 
Loading...
Searching...
No Matches
tstgmn.for
Go to the documentation of this file.
1C JJV changed name to ONECOV to avoid confusion with array COVAR
2C JJV this was also changed in the body of the function
3C REAL FUNCTION covar(x,y,n)
4 REAL function onecov(x,y,n)
5C .. Scalar Arguments ..
6 INTEGER n
7C ..
8C .. Array Arguments ..
9 REAL x(n),y(n)
10C ..
11C .. Local Scalars ..
12 REAL avx,avy,varx,vary,xmax,xmin
13 INTEGER i
14C ..
15C .. External Subroutines ..
16 EXTERNAL stat
17C ..
18C .. Intrinsic Functions ..
19 INTRINSIC real
20C ..
21C .. Executable Statements ..
22 CALL stat(x,n,avx,varx,xmin,xmax)
23 CALL stat(y,n,avy,vary,xmin,xmax)
24C covar = 0.0
25 onecov = 0.0
26 DO 10,i = 1,n
27C covar = covar + (x(i)-avx)* (y(i)-avy)
28 onecov = onecov + (x(i)-avx)* (y(i)-avy)
29 10 CONTINUE
30C covar = covar/real(n-1)
31 onecov = onecov/real(n-1)
32 RETURN
33
34 END
35
36C JJV Added argument LDXCOV (leading dimension of XCOVAR) to be
37C JJV consistent with the program TSTGMN, see comments below.
38C JJV This change necessitated changes in the declarations.
39C SUBROUTINE prcomp(p,mean,xcovar,answer)
40 SUBROUTINE prcomp(p,mean,xcovar,ldxcov,answer)
41
42C INTEGER p,maxp
43 INTEGER p,maxp,ldxcov
44 parameter(maxp=10)
45C REAL mean(p),xcovar(p,p),rcovar(maxp,maxp)
46 REAL mean(p),xcovar(ldxcov,p),rcovar(maxp,maxp)
47 REAL answer(1000,maxp)
48C JJV added ONECOV because of name change to function COVAR
49C 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
64C JJV changed COVAR to match new name
65C 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
75C JJV added LDCOV (leading dimension of COVAR) to be
76C JJV consistent with the program TSTGMN, see comments below.
77C JJV This change necessitated changes in the declarations.
78C SUBROUTINE setcov(p,var,corr,covar)
79 SUBROUTINE setcov(p,var,corr,covar,ldcov)
80C Set covariance matrix from variance and common correlation
81C .. Scalar Arguments ..
82 REAL corr
83C INTEGER p
84 INTEGER p,ldcov
85C ..
86C .. Array Arguments ..
87C REAL covar(p,p),var(p)
88 REAL covar(ldcov,p),var(p)
89C ..
90C .. Local Scalars ..
91 INTEGER i,j
92C ..
93C .. Intrinsic Functions ..
94 INTRINSIC sqrt
95C ..
96C .. 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)
112C .. Scalar Arguments ..
113 REAL av,var,xmax,xmin
114 INTEGER n
115C ..
116C .. Array Arguments ..
117 REAL x(n)
118C ..
119C .. Local Scalars ..
120 REAL sum
121 INTEGER i
122C ..
123C .. Intrinsic Functions ..
124 INTRINSIC real
125C ..
126C .. 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
146C Test Generation of Multivariate Normal Data
147C JJV SETGMN was: SUBROUTINE setgmn(meanv,covm,p,parm)
148C JJV is: SUBROUTINE setgmn(meanv,covm,ldcovm,p,parm)
149C JJV So the covariance matrices have been changed to 2-dim'l
150C JJV matrices, and the additional argument has been added to
151C JJV the subroutine call. Additional changes have been made
152C JJV to reflect this. (in declarations, the matrix copy routine,
153C JJV and in subroutine calls.)
154C .. Parameters ..
155 INTEGER maxp
156 parameter(maxp=10)
157 INTEGER maxobs
158 parameter(maxobs=1000)
159C JJV this parameter is no longer needed
160C INTEGER p2
161C PARAMETER (p2=maxp*maxp)
162C ..
163C .. Local Scalars ..
164 REAL corr
165 INTEGER i,iobs,is1,is2,j,p
166 CHARACTER phrase*100
167C ..
168C .. Local Arrays ..
169C REAL answer(1000,maxp),ccovar(p2),covar(p2),mean(maxp),param(500),
170C + 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)
173C ..
174C .. External Subroutines ..
175 EXTERNAL genmn,phrtsd,prcomp,setall,setcov,setgmn
176C ..
177C .. 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
197C 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)
203C DO 20,i = 1,p2
204C ccovar(i) = covar(i)
205C 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
211C
212C Generate Variables
213C
214C 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
222C CALL prcomp(p,mean,covar,answer)
223 CALL prcomp(p,mean,covar,maxp,answer)
224C
225C Print Comparison of Generated and Reconstructed Values
226C
227 GO TO 10
228
229 END
real
ColumnVector real(const ComplexColumnVector &a)
Definition dColVector.cc:137
genmn
subroutine genmn(parm, x, work)
Definition genmn.f:2
xmin
octave_int< T > xmin(const octave_int< T > &x, const octave_int< T > &y)
Definition oct-inttypes.h:1316
xmax
octave_int< T > xmax(const octave_int< T > &x, const octave_int< T > &y)
Definition oct-inttypes.h:1306
phrtsd
subroutine phrtsd(phrase, seed1, seed2)
Definition phrtsd.f:2
setall
subroutine setall(iseed1, iseed2)
Definition setall.f:2
setgmn
subroutine setgmn(meanv, covm, ldcovm, p, parm)
Definition setgmn.f:2
prcomp
subroutine prcomp(p, mean, xcovar, ldxcov, answer)
Definition tstgmn.for:41
setcov
subroutine setcov(p, var, corr, covar, ldcov)
Definition tstgmn.for:80
onecov
real function onecov(x, y, n)
Definition tstgmn.for:5
tstgmn
program tstgmn
Definition tstgmn.for:145
stat
subroutine stat(x, n, av, var, xmin, xmax)
Definition tstgmn.for:112