program DGMV
!_______________________________________________________________________
! A driver for the programs required for multivariate statisitics
! spring 2002.
!             Version 1/23/02            DRH
!_______________________________________________________________________
USE DFIMSLSS
USE IMSLf90 
implicit none
integer, parameter :: num=10000, p=3
integer :: i,j,k,n(5),l,cycles,h,g 
real :: z1(num),z2(num),z3(num),y(num,1),zy(num,1),zX(num,p)
real :: x1(num),x2(num),x3(num),avgx1,avgx2,avgx3,newy(num,1)
real :: varcovX(p,p),varx1,varx2,varx3,avgy,vary,cor,zXt(num,3)
real :: corX(p,p),Bt(p,1),B(p,1),X(num,p),corxy(p,1),Bn(p,1) 
real :: covX,Exp_Cov(p,p),avg_sum,sum,CovT(p)
data  n /10,50,100,500,1000/  !sample size vector


open (4,file='DGMV.txt')


!This program is designed to test how sensitive a covarience matrix
!is to sample size.  Since we can compute the expected values of a 
!covarience matrix when the variences of the variables put into are 
!known, we can calculate the difference between what we expect and 
!what we see when we test various sample sizes.  All the z variables 
!here are normal(0,1) and the x's are linear combinations of those.
!The expected covarience matrix is (remember covarience matrices are 
!the same on top and bottom):

!(var(x1)+var(x2)=1.0+1.0=2.0) , (var(x2)=1.0), (var(x1)=1.0) 
!				(var(x2)+var(x3)=1.0+1.0=2.0) , (var(x3)=1.0)
!								(var(x1)+var(x3)=1.0+1.0=2.0)

!This program will compute the summed absolute value of the difference
!between the observed and expected elements of the respective cov 
!matrices and average them for 100 trials at each sample size, 10,50
!100,500 and 100 and print the results to an output file.
do i=1,p
	do j=1,p
		if(i==j) then
			Exp_Cov(i,j)=2.0
		else
			Exp_Cov(i,j)=1.0
		end if
	end do
end do

!header for output file
write(4,'(a,a)') 'Sample Size  ','  Average sum absolute deviation from expected' 
cycles = 100 
do l=1,5  !cycle through each of the possible samle sizes
	do i=1,cycles !sample 100 times

	!RandomNumberGenerator(Typ,high,low,n,res)
	!for a uniform r.n. specify upper and lower bounds with high and low
	!for a normal, specify desired mean and std. dev. with high and low,
	!respectively.  *******BOTH MUST REAL******
		
		call RandomNumberGenerator('n',0.0,1.0,n(l),z1)
		call RandomNumberGenerator('n',0.0,1.0,n(l),z2)
		call RandomNumberGenerator('n',0.0,1.0,n(l),z3)
		do j=1,n(l) !create X vector, linear combo of z's
			!x1(j)=z1(i)+z2(i)
			X(j,1)=z1(j)+z2(j)
			X(j,2)=z2(j)+z3(j)
			X(j,3)=z1(j)+z3(j)
		end do
						
		!compute covarience matrix
		sum = 0.0
		do j=1,p
			do k=1,p
				call covarience(n(l),X(:,j),X(:,k),covX)
				sum = sum + abs(covX-Exp_Cov(j,k))
			end do
		end do
		avg_sum = avg_sum + sum
	end do
	avg_sum = avg_sum/real(cycles)
	write(4,'(i4,10x,F12.7)') n(l),avg_sum
	avg_sum =0.0
end do

end program DGMV