Fixed Model

 

자료

 

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 effects와 observations

 

프로그램은 변형된 LHS와 RHS를 구성하는 프로그램과 정렬된 LHS를 이용하여 해를 찾는 프로그램으로 구분

 

1) 변형된 LHS와 RHS를 구성하는 프로그램

 

PROGRAM Fixed_Model_Setup

! program name - fixed model setup
! programmer - Park Byoungho
! usage - fixed_model_setup < parameter_file
! purpose : read data file and write left-hand side and right-hand side
! Date - 2009.3.10.
! update - none

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, ALLOCATABLE :: fixed_data(:) ! fixed effects for each line
REAL :: value ! value of trait for each line

REAL, ALLOCATABLE :: xpy(:) ! right-hand side

INTEGER :: No_Of_Equations ! number of equation
INTEGER :: i, j ! do loop
REAL :: temp_value ! temporary storate for swap

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

! initialize xpy - right-hand side
No_Of_Equations = SUM(No_Of_Levels)
ALLOCATE(xpy(No_Of_Equations))
xpy = 0 ! initialize

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

IF (status == 0) THEN ! file open success

! open file for writing left-hand side
OPEN(UNIT = 2, FILE = 'lhs.dat', STATUS = 'REPLACE', ACTION = 'WRITE')

! read each line
DO
READ(UNIT = 1, FMT = *, IOSTAT = status) fixed_data, value
IF (status /= 0) EXIT ! end of file

! find location of equations
DO i = 2, No_Of_Effects
fixed_data(i) = fixed_data(i) + SUM(No_Of_Levels(1:i - 1))
END DO

! write location of non-zero elements :
DO i = 1, No_Of_Effects
! swap
temp_value = fixed_data(i)
fixed_data(i) = fixed_data(1)
fixed_data(1) = temp_value

WRITE(UNIT = 2, FMT = *) fixed_data, value
END DO

! add value to xpy
DO i = 1, No_Of_Effects
xpy(fixed_data(i)) = xpy(fixed_data(i)) + value
END DO

END DO ! end of reading data

CLOSE(1) ! close data file
CLOSE(2) ! close file for left-hand side

! open file for right hand side
OPEN(UNIT = 3, FILE = 'rhs.dat', STATUS = 'REPLACE', ACTION = 'WRITE')

! write xpy
DO i = 1, No_Of_Equations
WRITE(UNIT = 3, FMT = *) xpy(i)
END DO

CLOSE(3)! close file for right-hand side

! making parameter file to need to execute fixed_model_solve.exe
OPEN(UNIT = 4, FILE = 'fixed_model_solve.par', STATUS = 'REPLACE', ACTION = 'WRITE')
WRITE(UNIT = 4, FMT = *) No_OF_Effects
WRITE(UNIT = 4, FMT = *) No_Of_Equations
CLOSE(4)

ELSE ! data file open fail
WRITE (*,'(1X, A, A)') 'File open failed -- error message : ', error_msg
STOP
END IF

END PROGRAM Fixed_Model_Setup

 

위 프로그램의 결과 생긴 변형된 LHS

 

1 4 28.000000
4 1 28.000000
1 5 19.000000
5 1 19.000000
1 6 25.000000
6 1 25.000000
2 4 30.000000
4 2 30.000000
2 6 21.000000
6 2 21.000000
2 6 30.000000
6 2 30.000000
3 5 21.000000
5 3 21.000000
3 4 27.000000
4 3 27.000000
3 6 25.000000
6 3 25.000000
3 5 19.000000
5 3 19.000000

두번째 fixed effect를 equation number로 바꾸어 주고, 즉 (두번째 fixed effect + 첫번째 fixed effect의 레벨)하여

 

오리지날 LHS의 off-diagonal 의 위치를 찾는다.

 

위 프로그램에서 만들어진 RHS

 

72.000000
81.000000
92.000000
85.000000
59.000000
101.00000

 

2) sort

 

sort는 도스의 sort는 사용하지 말것. 각 라인을 문자열에 의하여 정렬하므로 소용없음

 

윈도우즈에 사용할 수 있는 유닉스 명령어 프로그램을 설치

 

http://blog.paran.com/bhpark70/30846064

 

정렬된 변형 LHS

 

1 4 28.000000
1 5 19.000000
1 6 25.000000
2 4 30.000000
2 6 21.000000
2 6 30.000000
3 4 27.000000
3 5 19.000000
3 5 21.000000
3 6 25.000000
4 1 28.000000
4 2 30.000000
4 3 27.000000
5 1 19.000000
5 3 19.000000
5 3 21.000000
6 1 25.000000
6 2 21.000000
6 2 30.000000
6 3 25.000000

첫 째 컬럼은 equation number

둘 째 컬럼은 해당 equation에서 0이 아닌 위치

 

 

 

3) 정렬된 변형 LHS와 RHS를 이용하여 해를 구하는 프로그램

 

PROGRAM fixed_model_solve

! program name - fixed_model_solve
! programmer - Park Byoungho
! usage - fixed_model_solve < fixed_model_solve.par
! purpose : read lhs_sorted.dat and rhs.dat which are made in fixed_model_setup.exe and are sorted
! find a solution using Gauss-Seidel iteration
! Date - 2009.3.15
! update - none

IMPLICIT NONE

! data dictionary
INTEGER :: no_of_effects ! number of effects
INTEGER :: pre_eq_no ! previous equation number
INTEGER :: cur_eq_no ! current equation number
INTEGER, ALLOCATABLE :: loc_in_eq(:) ! location of off-diagonal in the equation
REAL :: value ! observation
INTEGER :: level_of_eq ! diagonal element of x'x
INTEGER :: no_of_eq ! number of equations
REAL, ALLOCATABLE :: rhs(:) ! right-hand side, x'y
REAL :: rhs_element ! right-hand side one element
REAL, ALLOCATABLE :: solutions(:) ! solutions
REAL :: pre_sol ! previous solution
INTEGER :: iteration ! iteration
REAL :: epsilon ! sum of square between old and new solutions
INTEGER :: i ! loop
INTEGER :: status ! i/o status
CHARACTER(LEN = 40) :: error_msg ! error message
INTEGER, PARAMETER :: MAX_ITER = 10000 ! maximum number of iteration
REAL, PARAMETER :: criteria = 1.E-8 ! criteria for stopping

! get the number of effectss
WRITE (*,*) "Enter the number of effects : "
READ (*,*) no_of_effects

ALLOCATE(loc_in_eq(no_of_effects - 1))

! get the number of equations
WRITE (*,*) "Enter the number of equations : "
READ (*,*) no_of_eq

ALLOCATE(rhs(no_of_eq), solutions(no_of_eq))

! initialiaze solutions to zero
solutions = 0

!open right-hand side file
OPEN(UNIT = 11, FILE = "rhs.dat", STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg)

IF (status == 0) THEN ! file open success

! read and store right-hand side
READ(UNIT = 11,FMT = *) rhs

ELSE ! file open fail
WRITE (*,'(1X, A, A)') 'File open failed -- error message : ', error_msg
STOP
END IF

close(11)

!open sorted left-hand side file
OPEN(UNIT = 12, FILE = "sorted_lhs.dat", STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg)

IF (status == 0) THEN ! file open success

DO iteration = 1, MAX_ITER

pre_eq_no = 0
level_of_eq = 0
epsilon = 0

DO
! read each line
READ(UNIT = 12, FMT = *, IOSTAT = status) cur_eq_no, loc_in_eq, value

IF (status == 0) THEN ! read data normally

IF (cur_eq_no == pre_eq_no) THEN ! same equation number
rhs_element = rhs_element - sum(solutions(loc_in_eq))
level_of_eq = level_of_eq + 1
pre_eq_no = cur_eq_no
ELSE ! new equation number
IF (cur_eq_no > 1) THEN
pre_sol = solutions(pre_eq_no) ! store old solution
solutions(pre_eq_no) = rhs_element / level_of_eq ! get a new solution
epsilon = epsilon + (pre_sol - solutions(pre_eq_no)) ** 2 ! calculate sum of squares of difference between old solution and new solution
!WRITE(*,*) pre_eq_no, level_of_eq, solutions(pre_eq_no)
END IF

pre_eq_no = cur_eq_no
rhs_element = rhs(cur_eq_no) - sum(solutions(loc_in_eq))
level_of_eq = 1
END IF

ELSE ! end of file
pre_sol = solutions(pre_eq_no) ! store old solution
solutions(pre_eq_no) = rhs_element / level_of_eq ! get a new solution
epsilon = epsilon + (pre_sol - solutions(pre_eq_no)) ** 2 ! calculate sum of squares of difference between old solution and new solutions
!WRITE(*,*) cur_eq_no, level_of_eq, solutions(pre_eq_no)
EXIT
END IF

END DO ! one iteration ends

! write iteration and epsilon
epsilon = epsilon / no_of_eq
WRITE(*,*) 'iteration = ', iteration , ', epsilon = ', epsilon
IF (epsilon < criteria) THEN
EXIT
ELSE
REWIND(12)
END IF

END DO

CLOSE(12)

! open file for writing solutions
OPEN(UNIT = 13, FILE = 'sol.dat', STATUS = 'REPLACE', ACTION = 'WRITE', IOSTAT = status, IOMSG = error_msg)
DO i = 1, no_of_eq
WRITE(UNIT = 13, FMT = *) solutions(i)
END DO
CLOSE(13)

ELSE ! sorted left-hand side file open failed
WRITE (*,'(1X, A, A)') 'File open failed -- error message : ', error_msg
STOP
END IF

 

END PROGRAM fixed_model_solve

 

위 프로그램을 이용하여 구한 solution

 

24.054985
25.405043
24.154963
3.7950027
-4.4549699
0.49499130