GNU Octave  4.0.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
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
setcov
subroutine setcov(p, var, corr, covar, ldcov)
Definition: tstgmn.for:79
phrtsd
subroutine phrtsd(phrase, seed1, seed2)
Definition: phrtsd.f:1
onecov
real function onecov(x, y, n)
Definition: tstgmn.for:4
tstgmn
program tstgmn
Definition: tstgmn.for:145
prcomp
subroutine prcomp(p, mean, xcovar, ldxcov, answer)
Definition: tstgmn.for:40
setgmn
subroutine setgmn(meanv, covm, ldcovm, p, parm)
Definition: setgmn.f:1
setall
subroutine setall(iseed1, iseed2)
Definition: setall.f:1
stat
subroutine stat(x, n, av, var, xmin, xmax)
Definition: tstgmn.for:111
genmn
subroutine genmn(parm, x, work)
Definition: genmn.f:1
real
ColumnVector real(const ComplexColumnVector &a)
Definition: dColVector.cc:156
octave_value::sqrt
octave_value sqrt(void) const
Definition: ov.h:1200
Generated on Thu Jun 4 2015 23:30:24 for GNU Octave by   doxygen 1.8.8