home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / libcruft / ranlib / tstgmn.for < prev    next >
Text File  |  1996-09-28  |  5KB  |  186 lines
  1.       REAL FUNCTION covar(x,y,n)
  2. C     .. Scalar Arguments ..
  3.       INTEGER n
  4. C     ..
  5. C     .. Array Arguments ..
  6.       REAL x(n),y(n)
  7. C     ..
  8. C     .. Local Scalars ..
  9.       REAL avx,avy,varx,vary,xmax,xmin
  10.       INTEGER i
  11. C     ..
  12. C     .. External Subroutines ..
  13.       EXTERNAL stat
  14. C     ..
  15. C     .. Intrinsic Functions ..
  16.       INTRINSIC real
  17. C     ..
  18. C     .. Executable Statements ..
  19.       CALL stat(x,n,avx,varx,xmin,xmax)
  20.       CALL stat(y,n,avy,vary,xmin,xmax)
  21.       covar = 0.0
  22.       DO 10,i = 1,n
  23.           covar = covar + (x(i)-avx)* (y(i)-avy)
  24.    10 CONTINUE
  25.       covar = covar/real(n-1)
  26.       RETURN
  27.  
  28.       END
  29.       SUBROUTINE prcomp(p,mean,xcovar,answer)
  30.  
  31.       INTEGER p,maxp
  32.       PARAMETER (maxp=10)
  33.       REAL mean(p),xcovar(p,p),rcovar(maxp,maxp)
  34.       REAL answer(1000,maxp)
  35.       REAL rmean(maxp),rvar(maxp)
  36.       INTEGER maxobs
  37.       PARAMETER (maxobs=1000)
  38.  
  39.       DO 10,i = 1,p
  40.           CALL stat(answer(1,i),maxobs,rmean(i),rvar(i),dum1,dum2)
  41.           WRITE (*,*) ' Variable Number',i
  42.           WRITE (*,*) ' Mean ',mean(i),' Generated ',rmean(i)
  43.           WRITE (*,*) ' Variance ',xcovar(i,i),' Generated',rvar(i)
  44.    10 CONTINUE
  45.       WRITE (*,*) '                   Covariances'
  46.       DO 30,i = 1,p
  47.           DO 20,j = 1,i - 1
  48.               WRITE (*,*) ' I = ',i,' J = ',j
  49.               rcovar(i,j) = covar(answer(1,i),answer(1,j),maxobs)
  50.               WRITE (*,*) ' Covariance ',xcovar(i,j),' Generated ',
  51.      +          rcovar(i,j)
  52.    20     CONTINUE
  53.    30 CONTINUE
  54.       RETURN
  55.  
  56.       END
  57.       SUBROUTINE setcov(p,var,corr,covar)
  58. C     Set covariance matrix from variance and common correlation
  59. C     .. Scalar Arguments ..
  60.       REAL corr
  61.       INTEGER p
  62. C     ..
  63. C     .. Array Arguments ..
  64.       REAL covar(p,p),var(p)
  65. C     ..
  66. C     .. Local Scalars ..
  67.       INTEGER i,j
  68. C     ..
  69. C     .. Intrinsic Functions ..
  70.       INTRINSIC sqrt
  71. C     ..
  72. C     .. Executable Statements ..
  73.       DO 40,i = 1,p
  74.           DO 30,j = 1,p
  75.               IF (.NOT. (i.EQ.j)) GO TO 10
  76.               covar(i,j) = var(i)
  77.               GO TO 20
  78.  
  79.    10         covar(i,j) = corr*sqrt(var(i)*var(j))
  80.    20         CONTINUE
  81.    30     CONTINUE
  82.    40 CONTINUE
  83.       RETURN
  84.  
  85.       END
  86.       SUBROUTINE stat(x,n,av,var,xmin,xmax)
  87. C     .. Scalar Arguments ..
  88.       REAL av,var,xmax,xmin
  89.       INTEGER n
  90. C     ..
  91. C     .. Array Arguments ..
  92.       REAL x(n)
  93. C     ..
  94. C     .. Local Scalars ..
  95.       REAL sum
  96.       INTEGER i
  97. C     ..
  98. C     .. Intrinsic Functions ..
  99.       INTRINSIC real
  100. C     ..
  101. C     .. Executable Statements ..
  102.       xmin = x(1)
  103.       xmax = x(1)
  104.       sum = 0.0
  105.       DO 10,i = 1,n
  106.           sum = sum + x(i)
  107.           IF (x(i).LT.xmin) xmin = x(i)
  108.           IF (x(i).GT.xmax) xmax = x(i)
  109.    10 CONTINUE
  110.       av = sum/real(n)
  111.       sum = 0.0
  112.       DO 20,i = 1,n
  113.           sum = sum + (x(i)-av)**2
  114.    20 CONTINUE
  115.       var = sum/real(n-1)
  116.       RETURN
  117.  
  118.       END
  119.       PROGRAM tstgmn
  120. C     Test Generation of Multivariate Normal Data
  121. C     .. Parameters ..
  122.       INTEGER maxp
  123.       PARAMETER (maxp=10)
  124.       INTEGER maxobs
  125.       PARAMETER (maxobs=1000)
  126.       INTEGER p2
  127.       PARAMETER (p2=maxp*maxp)
  128. C     ..
  129. C     .. Local Scalars ..
  130.       REAL corr
  131.       INTEGER i,iobs,is1,is2,j,p
  132.       CHARACTER phrase*100
  133. C     ..
  134. C     .. Local Arrays ..
  135.       REAL answer(1000,maxp),ccovar(p2),covar(p2),mean(maxp),param(500),
  136.      +     temp(maxp),var(maxp),work(maxp)
  137. C     ..
  138. C     .. External Subroutines ..
  139.       EXTERNAL genmn,phrtsd,prcomp,setall,setcov,setgmn
  140. C     ..
  141. C     .. Executable Statements ..
  142.       WRITE (*,9000)
  143.  
  144.  9000 FORMAT (
  145.      +     ' Tests Multivariate Normal Generator for Up to 10 Variables'
  146.      +       /
  147.      +  ' User inputs means, variances, one correlation that is applied'
  148.      +       /'     to all pairs of variables'/
  149.      +       ' 1000 multivariate normal deviates are generated'/
  150.      +     ' Means, variances and covariances are calculated for these.'
  151.      +       )
  152.  
  153.    10 WRITE (*,*) 'Enter number of variables for normal generator'
  154.       READ (*,*) p
  155.       WRITE (*,*) 'Enter mean vector of length ',p
  156.       READ (*,*) (mean(i),i=1,p)
  157.       WRITE (*,*) 'Enter variance vector of length ',p
  158.       READ (*,*) (var(i),i=1,p)
  159.       WRITE (*,*) 'Enter correlation of all variables'
  160.       READ (*,*) corr
  161.       CALL setcov(p,var,corr,covar)
  162.       WRITE (*,*) ' Enter phrase to initialize rn generator'
  163.       READ (*,'(a)') phrase
  164.       CALL phrtsd(phrase,is1,is2)
  165.       CALL setall(is1,is2)
  166.       DO 20,i = 1,p2
  167.           ccovar(i) = covar(i)
  168.    20 CONTINUE
  169. C
  170. C     Generate Variables
  171. C
  172.       CALL setgmn(mean,ccovar,p,param)
  173.       DO 40,iobs = 1,maxobs
  174.           CALL genmn(param,work,temp)
  175.           DO 30,j = 1,p
  176.               answer(iobs,j) = work(j)
  177.    30     CONTINUE
  178.    40 CONTINUE
  179.       CALL prcomp(p,mean,covar,answer)
  180. C
  181. C     Print Comparison of Generated and Reconstructed Values
  182. C
  183.       GO TO 10
  184.  
  185.       END
  186.