113 lines
3.7 KiB
FortranFixed
113 lines
3.7 KiB
FortranFixed
! program main
|
|
! implicit none
|
|
! integer nsamp,mdim,mpc
|
|
! double precision sample(100,100),princomp(100,100),
|
|
! &transdata(100,100),x(100),eigenvector(100,100),eigenvalue(100)
|
|
! integer i,j,k
|
|
! open(unit=1,file='Table8.3.txt')
|
|
! read(1,*)
|
|
! nsamp=0
|
|
!10 read(1,*,end=100)i,(x(j),j=1,6)
|
|
! if(i.le.1)goto 10
|
|
! nsamp=nsamp+1
|
|
! do j=1,6
|
|
! sample(nsamp,j)=x(j)
|
|
! enddo
|
|
! goto 10
|
|
!100 close(1)
|
|
! mdim=6
|
|
! mpc=2
|
|
! call princompana(nsamp,mdim,sample(1:nsamp,1:mdim),mpc,
|
|
! &princomp(1:nsamp,1:mpc),transdata(1:nsamp,1:mdim),
|
|
! &eigenvector(1:mdim,1:mdim),eigenvalue)
|
|
! do i=1,mpc
|
|
! do j=1,nsamp
|
|
! write(*,*)j,princomp(j,i)
|
|
! enddo
|
|
! enddo
|
|
! do i=1,mdim
|
|
! do j=1,nsamp
|
|
! write(*,*)j,transdata(j,i)
|
|
! enddo
|
|
! enddo
|
|
! end
|
|
|
|
subroutine princompana(nsamp,mdim,sample,mpc,princomp,transdata,
|
|
&eigenvector,eigenvalue)
|
|
implicit none
|
|
!-----------Inputs----------------------------------------
|
|
!mpc is the number of principal components to keep
|
|
integer nsamp,mdim,mpc
|
|
double precision sample(nsamp,mdim)
|
|
!-----------Outputs---------------------------------------
|
|
!princomp is the projection of a sample on the principal axes
|
|
!transdata is the data of the orginal sample filtered with mpc principal components
|
|
double precision eigenvalue(mdim),eigenvector(mdim,mdim),
|
|
&sampmean(mdim),princomp(nsamp,mpc),transdata(nsamp,mdim),
|
|
&sampadj(nsamp,mdim)
|
|
!---------------------------------------------------------
|
|
integer i,j,k
|
|
call geteigen(nsamp,mdim,sample(1:nsamp,1:mdim),eigenvalue,
|
|
&eigenvector(1:mdim,1:mdim),sampmean,sampadj(1:nsamp,1:mdim))
|
|
do i=1,mpc
|
|
do j=1,nsamp
|
|
princomp(j,i)=0.0d0
|
|
do k=1,mdim
|
|
princomp(j,i)=princomp(j,i)+eigenvector(k,i)*sampadj(j,k)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
do j=1,mdim
|
|
do i=1,nsamp
|
|
transdata(i,j)=sampmean(j)
|
|
do k=1,mpc
|
|
transdata(i,j)=transdata(i,j)+eigenvector(j,k)*princomp(i,k)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
return
|
|
end
|
|
|
|
subroutine geteigen(nsamp,mdim,sample,eigenvalue,eigenvector,
|
|
&sampmean,sampadj)
|
|
integer nsamp,mdim
|
|
double precision sample(nsamp,mdim),eigenvalue(mdim),
|
|
&eigenvector(mdim,mdim),sampmean(mdim),sampadj(nsamp,mdim)
|
|
!Each column is an eigenvector. The first column corresponds to the largest eigenvalue and the last column corresponds
|
|
!to the smallest eigenvalue
|
|
call covariancematrix(nsamp,mdim,sample(1:nsamp,1:mdim),
|
|
&eigenvector(1:mdim,1:mdim),sampmean,sampadj(1:nsamp,1:mdim))
|
|
call eigen_sym_up(mdim,eigenvector(1:mdim,1:mdim),eigenvalue)
|
|
return
|
|
end
|
|
|
|
subroutine covariancematrix(nsamp,mdim,sample,covarmatrix,
|
|
&sampmean,sampadj)
|
|
implicit none
|
|
!covarmatrix is an upper trangle
|
|
integer nsamp,mdim
|
|
double precision sample(nsamp,mdim),covarmatrix(mdim,mdim),
|
|
&sampmean(mdim),sampadj(nsamp,mdim)
|
|
integer i,j,k
|
|
do j=1,mdim
|
|
sampmean(j)=0.0d0
|
|
do i=1,nsamp
|
|
sampmean(j)=sampmean(j)+sample(i,j)
|
|
enddo
|
|
sampmean(j)=sampmean(j)/dble(nsamp)
|
|
do i=1,nsamp
|
|
sampadj(i,j)=sample(i,j)-sampmean(j)
|
|
enddo
|
|
enddo
|
|
do i=1,mdim
|
|
do j=i,mdim
|
|
covarmatrix(i,j)=0.0d0
|
|
do k=1,nsamp
|
|
covarmatrix(i,j)=covarmatrix(i,j)+sampadj(k,i)*sampadj(k,j)
|
|
enddo
|
|
covarmatrix(i,j)=covarmatrix(i,j)/dble(nsamp-1)
|
|
enddo
|
|
enddo
|
|
return
|
|
end
|