Subroutine average(n,vec,avg)
implicit none
integer :: n,i
real :: vec(n),avg
avg=0.0
do i=1,n
	avg = avg + vec(i)
end do
avg = avg/real(n)
end subroutine average
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subroutine varience(n,vec,var)
implicit none
integer :: n,i
real :: vec(n),avg,var
call average(n,vec,avg)
var=0.0
do i=1,n
	var = var + ((vec(i)-avg)**2.0)
end do
var = var/real(n-1)
end subroutine varience
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subroutine variance(n,vec,var)
implicit none
integer :: n,i
real :: vec(n),avg,var
call average(n,vec,avg)
var=0.0
do i=1,n
	var = var + ((vec(i)-avg)**2.0)
end do
var = var/real(n-1)
end subroutine variance
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subroutine covarience(n,x,y,cov)
implicit none
integer :: n,i
real :: x(n),y(n),xy(n),cov,avgx,avgy,avgxy
avgx=0.0
avgy=0.0
avgxy=0.0
call average(n,x,avgx)
call average(n,y,avgy)
do i=1,n
	xy(i) = x(i) * y(i)
enddo
call average(n,xy,avgxy)
cov=avgxy-avgx*avgy
end subroutine covarience
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subroutine covariance(n,x,y,cov)
implicit none
integer :: n,i
real :: x(n),y(n),xy(n),cov,avgx,avgy,avgxy
avgx=0.0
avgy=0.0
avgxy=0.0
call average(n,x,avgx)
call average(n,y,avgy)
do i=1,n
	xy(i) = x(i) * y(i)
enddo
call average(n,xy,avgxy)
cov=avgxy-avgx*avgy
end subroutine covariance
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subroutine correlation(n,x,y,cor)
implicit none
integer :: n,i
real :: x(n),y(n),cov,cor,vary,varx
cor=0.0
cov=0.0
call covarience(n,x,y,cov)
call varience(n,x,varx)
call varience(n,y,vary)
cor = cov/sqrt(varx*vary)
end subroutine correlation
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
integer function fac(n)
implicit none
integer :: n,i
fac=1.0
do i=n,1,i-1
	fac=fac*i
end do
end function fac
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
integer function choose(n,r)
implicit none
integer :: n,r
integer, external :: fac
choose = fac(n)/(fac(r)*fac(n-r))
end function choose
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%