Fixed Model 2

자료로부터 X, y X'X, X'y, X'X의 일반화 역행렬, solution을 구하는 프로그램

 

 

자료

 

1 1 28
1 2 19
1 3 25
2 1 30
2 3 21
2 3 30
3 2 21
3 1 27
3 3 25
3 2 19
위 자료를 fixed_model.dat로 저장

 

 

프로그램 소스

 

PROGRAM Fixed_Model
! program name - fixed_model
! programmer - Park Byoungho
! usage - fixed_model < parameter_file
! purpose : read data file, construct left-hand side and right-hand side and find solution
! Date - 2009.4.7.
! update - none USE gi

IMPLICIT NONE

! data dictionary
CHARACTER(LEN = 40) :: filename ! data file name
INTEGER :: No_Of_Effects ! number of effects
INTEGER, ALLOCATABLE :: No_Of_Levels(:) ! number of levels of each fixed effect
INTEGER :: no_of_obs ! number of observations

INTEGER, ALLOCATABLE :: fixed_data(:) ! fixed effects for each line REAL, ALLOCATABLE :: x(:,:), y(:) ! read and store data from the file
REAL, ALLOCATABLE :: xp(:,:), xpx(:,:), xpy(:), gxpx(:,:), sol(:) ! transpose of x, left-hand side, right-hand side

INTEGER :: No_Of_Equations ! number of equation
INTEGER :: i, j ! do loop

INTEGER :: status ! I/O status
CHARACTER(LEN = 40) :: error_msg ! error message

! get data file name
WRITE (*,*) "Enter data file name : "
READ (*,*) filename ! get number of fixed effects
WRITE (*,*) "Enter the number of fixed effects : "
READ (*,*) No_Of_Effects ! allocate dimensions of no_of_levels and fixed_data
ALLOCATE(No_Of_Levels(No_Of_Effects), fixed_data(No_Of_Effects))

WRITE (*,*) "Enter the number of levels for each fixed effect : "
READ (*,*) No_Of_Levels ! get number of observations
WRITE (*,*) "Enter the number of observations : "
READ (*,*) no_of_obs ! initialize xpy - right-hand side
no_of_equations = SUM(no_of_levels)
ALLOCATE(x(no_of_obs,no_of_equations),&
y(no_of_obs),&
xp(no_of_equations,no_of_obs),&
xpx(no_of_equations,no_of_equations), &
gxpx(no_of_equations,no_of_equations), &
xpy(no_of_equations), &
sol(no_of_equations) &
) ! initialize
x = 0
y = 0
xp = 0
xpx = 0
xpy = 0
gxpx = 0

! open data file
OPEN(UNIT = 11, FILE = filename, STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg)

IF (status /= 0) THEN ! file open failed
WRITE (*,'(1X, A, A)') 'File open failed -- error message : ', error_msg
STOP
END IF DO i = 1, no_of_obs ! read each line
READ(UNIT = 11, FMT = *, IOSTAT = status) fixed_data, y(i)
IF (status /= 0) EXIT ! end of file ! find location of equations
DO j = 2, No_Of_Effects
fixed_data(j) = fixed_data(j) + SUM(no_of_levels(1:j - 1))
END DO

x(i,fixed_data) = 1 END DO xp = TRANSPOSE(x)
xpx = MATMUL(xp, x)
xpy = MATMUL(xp, y) CALL geninv(xpx, gxpx)

sol = MATMUL(gxpx, xpy)

WRITE(*,*) 'X ='
WRITE(*,'(6(F3.0, 1X))') ((x(i,j), j = 1, no_of_equations), i = 1, no_of_obs)

WRITE(*,*) 'y ='
WRITE(*,'(F6.1)') (y(i), i = 1, no_of_obs)

WRITE(*,*) 'X'' = '
WRITE(*,'(10(F3.0, 1X))') ((xp(i,j),j = 1, no_of_obs), i = 1, no_of_equations)

WRITE(*,*) 'X''X ='
WRITE(*,'(6(F3.0, 1X))') ((xpx(i,j), j = 1, no_of_equations), i = 1, no_of_equations)

WRITE(*,*) 'X''y ='
WRITE(*,'(F6.1)') (xpy(i),i = 1, no_of_equations) WRITE(*,*) 'ginverse of X''X ='
WRITE(*,'(6(F6.3, 1X))') ((gxpx(i,j), j = 1, no_of_equations), i = 1, no_of_equations)

WRITE(*,*) 'solutions ='
WRITE(*,'(F6.1)') (sol(i),i = 1, no_of_equations)

END PROGRAM Fixed_Model

 

위 프로그램을 fixed_model.f95로 저장

 

프로그램을 편하게 실행하기 위하여fixed_model.dat
2
3 3
10
위 자료를 fixed_model.par로 저장

 

일반화 역행렬(generalized inverse)을 구하는 모듈

 

MODULE gi CONTAINS

SUBROUTINE geninv(x,gix) REAL, DIMENSION(:,:), INTENT(IN) :: x ! square matrix to inverse
REAL, DIMENSION(:,:), INTENT(OUT) :: gix ! square matirx after inverse
INTEGER :: i, j, no_of_eq, size_x, mr
REAL(KIND = 8), ALLOCATABLE :: upperd(:)
INTEGER, DIMENSION(99) :: iw
REAL(KIND = 8) :: z no_of_eq = SIZE(x,1)

ALLOCATE(upperd((no_of_eq*(no_of_eq +1)/2)))

upperd = 0

DO i = 1, no_of_eq
DO j = 1, no_of_eq
upperd(ihmssf(i,j,no_of_eq))=x(i,j)
END DO
END DO CALL DJNVHF(upperd,no_of_eq,iw,z,mr) DO i = 1, no_of_eq
DO j = 1, no_of_eq
gix(i,j) = upperd(ihmssf(i,j,no_of_eq))
END DO
END DO END SUBROUTINE geninv SUBROUTINE DJNVHF(A,N,IS,Z,IR)
!
! MATRIX INVERSION SUBROUTINE
!
! 'A' IS A HALF-STORED MATRIX TO BE INVERTED. THE INVERTED
! RESULT IS RETURNED IN 'A'.
! 'N' IS THE ORDER OF THE MATRIX TO BE INVERTED.
! 'IS' IS A WORK VECTOR.
! 'Z' IS SET TO ZERO ON CALLS TO THIS SUBROUTINE FROM THE BEEF
! SIRE MONITORING SYSTEM.
! 'IR' RETURNS THE RANK OF THE MATRIX.
!
! THIS SUBROUTINE CALLS THE SUBROUTINE DREARN.
!
DIMENSION A(1),IS(1)
REAL*8 A,BIG,Z,RECIP,R
201 NP=(N*(N+1))/2
IF(N-1) 21,23,96
21 WRITE(6,22)
22 FORMAT(' N.LT.1')
23 IF(A(1)) 25,24,25
24 IR=0
26 RETURN
25 IR=1
A(1)=1.D0/A(1)
RETURN
96 DO 50 L=1,N
BIG=0.D1
DO 2 I=L,N
II=-(I*(I-3))/2+N*(I-1)
IF(DABS(A(II))-DABS(BIG)) 2,2,1
1 IS(L)=I
BIG=A(II)
2 CONTINUE
IF(DABS(BIG)-Z) 63,63,62
63 IR=L-1
IF(IR) 26,26,98
98 DO 64 I=1,L
LI=-(I*(I-1))/2+N*(I-1)
DO 64 J=L,N
LIJ=LI+J
64 A(LIJ)=-0.D1
IF(L-N) 66,69,67
67 STOP 67
66 LP1=L+1
DO 68 I=LP1,N
LI=-(I*(I-1))/2+N*(I-1)
DO 68 J=I,N
LIJ=LI+J
68 A(LIJ)=-0.D1
69 IF(IR-2) 56,27,27
56 CALL DREARN(A,N,IS,1)
DO 57 I=1,NP
57 A(I)=-A(I)
RETURN
27 LP1=L+1
DO 65 I=2,L
J=LP1-I
65 CALL DREARN(A,N,IS,J)
DO 101 I=1,NP
101 A(I)=-A(I)
RETURN
62 CALL DREARN(A,N,IS,L)
28 K3=N*(L-1)-(L*(L-3))/2
RECIP=1./A(K3)
A(K3)=-RECIP
DO 50 I=1,N
IF(I-L) 6,50,7
6 K11=N*(I-1)-(I*(I-3)/2)
K1=K11+L-I
GO TO 8
7 K11=N*(I-1)-(I*(I-3))/2
K1=K3+I-L
8 R=A(K1)*RECIP
301 DO 12 J=I,N
K4=K11+J-I
IF(J-L) 10,12,11
10 K5=N*(J-1)-(J*(J-3)/2)+L-J
GO TO 14
11 K5=K3+J-L
14 A(K4)=A(K4)-R*A(K5)
12 CONTINUE
A(K1)=R
50 CONTINUE
NP=(N*(N+1))/2
DO 100 I=1,NP
100 A(I)=-A(I)
NP1=N+1
DO 61 I=2,N
L=NP1-I
61 CALL DREARN(A,N,IS,L)
IR=N
RETURN
END SUBROUTINE

SUBROUTINE DREARN(A,N,IS,L)
!
! THIS IS A SUPPORT SUBROUTINE CALLED BY THE MATRIX INVERSION
! SUBROUTINE DJNVHF.
!
DIMENSION A(1),IS(1)
DOUBLE PRECISION A,SAVE
ISL=IS(L)
IF(ISL-L) 3,28,4
3 STOP 3
4 LM1=L-1
IF(LM1) 22,22,5
5 DO 21 I=1,LM1
IL=-(I*(I-1))/2+N*(I-1)+L
IISL=IL-L+ISL
SAVE=A(IL)
A(IL)=A(IISL)
21 A(IISL)=SAVE
22 LP1=L+1
ISLM1=ISL-1
IF(LP1-ISLM1) 23,23,25
23 DO 24 I=LP1,ISLM1
LI=-(L*(L-1))/2+N*(L-1)+I
IISL=-(I*(I-1))/2+N*(I-1)+ISL
SAVE=A(LI)
A(LI)=A(IISL)
24 A(IISL)=SAVE
25 ISLP1=ISL+1
IF(ISLP1-N) 26,26,38
26 DO 27 I=ISLP1,N
ISLI=-(ISL*ISLM1)/2+N*ISLM1+I
LI=-(L*LM1)/2+N*LM1+I
SAVE=A(ISLI)
A(ISLI)=A(LI)
27 A(LI)=SAVE
38 LL=-(L*(L-3))/2+N*(L-1)
ISLISL=-(ISL*(ISL-3))/2+ISLM1*N
SAVE=A(LL)
A(LL)=A(ISLISL)
A(ISLISL)=SAVE
28 RETURN
END SUBROUTINE FUNCTION IHMSSF(IR,IC,N)

JR=IR
JC=IC
IF(IR-IC)40,40,41
41 JC=IR
JR=IC
40 NR=JR-1
KK=(JR*NR)/2
KK=KK+NR*(N-NR)
KK=KK+JC-NR
IHMSSF=KK
RETURN
END FUNCTION

 

END MODULE gi

 

위 소스를 inverse.f95로 저장

 

 

프로그램 컴파일 및 실행