program data_by_cover_matrix
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
! This program is designed to generate random numbers that conform to
!a pre-specified covariance matrix.  I will use the tired old mortality
!dataset as a template for the covariance matrix.
!			version 2/8/02          DRH
!
!here are the steps:
!1) Specify a covariance matrix S
!2) do an eigan decomposition on it
!3) generate a sample of n realizations of xhat.  These are generated as
!  follows; generate random numbers with variances equal to l, the vector 
!  of eigan values that come from S.  Each number is an observation, each
!  new variance corresponds to a new variable.  Experiment with sample
!  size.  Now take these and form a new matrix Y.  Xhat'=Uy' where U is the
!  matrix of eigan vectors from S
!4) Calculate a new covariance matrix Sx
!		How close is Sx to S?
!5) Do a PCA on Sx
!		How close the vector of eigan values from this, lx to l?
!		How close is Ux to U?
!6) Compute the PC scores, Yhat=U'Xx
!		How close is Yhat to Y?
!________________________________________________________________________
USE DFIMSLSS
USE IMSLf90 
implicit none

integer, parameter :: rows=59  !max size of data arrays
integer, parameter :: p=7  !look at your data file and check the # of predictors
integer, parameter :: ss=1000 !sample size of generated data
!YOU NEED A BIG SAMPLE SIZE GENERATE DATA WITH A SPECIFIED COVARIANCE
!USING THIS METHOD.

integer, external :: choose

integer :: i,j,k,g,n,q,o,ioerror,trial(p),place,col
integer :: incl_vars(p)

real ::  X(rows,p),y(rows,1), yhat(ss,1), Xt(p,rows) !
real :: yavg,Xavg,covar,var
real :: L(p,1),U(p,p),Uinv(p,p),Ut(p,p),Sx(p,p),UhatT(p,p)
real :: Xnew(ss,p),Xhat(ss,p),Sxhat(p,p),Lhat(p,1),Uhat(p,p)
real :: XhatT(p,ss),PC(p,ss),PCT(ss,p),PCold(rows,p),XnewT(p,ss)




open (4,file='data_gen.txt') !output file
open (20,file='Mortality_dataset.txt',status='old') !datafile

!read in data
n=0
ioerror=0
do i=1,rows
	n=n+1	
	read (20,*,iostat=ioerror) y(i,1),(X(i,j),j=1,p)
	if(ioerror==-1) then
		n=n-1
		write(6,'(a,i3)') 'end of file reached, n= ',n
		exit
	else if (ioerror/=0) then
		write(6,'(a,i3)') 'read error!, line ',i
	endif
	!write(4,'(<p>(F9.2,1x))') (X(i,j),j=1,p)
enddo
!**********************************************************************
!Create and write covariance matrix to output file
write(4,'(a)') 'covariance matrix of x = '
do i=1,p
	do j=1,p
		call covarience(n,X(:,i),X(:,j),covar)
		Sx(i,j) = covar
	end do
	write(4,'(<p>(f10.2,1x))') Sx(i,:)
end do

!We need to do an eigan decomposition of Sx now.
!this routine is for the eigan decomposition of a symmetric matrix
call EVCSF(p,Sx,p,L,U,p)
!L is the vector of eigan values

!generate new random numbers with specified variance = L
do j=1, ss
	do i=1,p
		var=0.0
		call RandomNumberGenerator('n',0.0,L(i,1),1,Xnew(j,i))		
	end do
end do
	!now we have to convert to new (generated) x values
XnewT = .t. Xnew  !not trying to throw this in a single line calculation
!saves RAM.
XhatT = U .x. XnewT
Xhat = .t. XhatT


!make a covar matrix to test our generated sample
write(4,'(/,a)') "covariance matrix of generated x's = "
do i=1,p
	do j=1,p
		call covarience(ss,Xhat(:,i),Xhat(:,j),covar)
		Sxhat(i,j) = covar
	end do
	write(4,'(<p>(f10.2,1x))') Sxhat(i,:)
end do
!look at the difference between the specified S and the S from
!the sample we generated to = it.
write(4,'(/,a)') "matrix of difference between them = "
do i=1,p
		write(4,'(<p>(f10.2,1x))') ((Sxhat(i,j)-Sx(i,j)),j=1,p)
end do
!Now perform PCA on the new covariance matrix
call EVCSF(p,Sxhat,p,Lhat,Uhat,p)

!now check these results
write(4,'(/,a)') 'L old           L of sample'
do i=1,p
	write(4,'(2(f15.4,2x))') L(i,1), Lhat(i,1)
end do
write(4,'(/,a)') "eigan vectors of old x's = "
do i=1,p
	write(4,'(<p>(f10.2,1x))') U(i,:)
end do
write(4,'(/,a)') "eigan vectors of generated x's = "
do i=1,p
	write(4,'(<p>(f10.2,1x))') Uhat(i,:)
end do

!Compute the PCA scores
XhatT = .t. Xhat
UhatT = .t. Uhat

PC =   UhatT .x. XhatT
PCT = .t. PC
!now compare this to Xnew above

!DON'T WRITE THESE LINES UNLESS YOU WANT A HUGE OUTPUT FILE.  THEY MATCH
!I CHECKED (THEY MATCH AT LEAST TO WITHIN SAMPLING ERROR).

!write(4,'(/,a)') "the new PC scores are "
do i=1,ss
!	write(4,'(<p>(f10.3,1x))') (PCT(i,j),j=1,p)
end do
!write(4,'(/,a)') "the old generated variables are "
do i=1,ss
!	write(4,'(<p>(f10.3,1x))') (Xnew(i,j),j=1,p)
end do


	stop 'normal completion'
end program data_by_cover_matrix