Fixed Regression Model - Fortran 프로그램

 

Fixed Regression Model

 

자료

 

1 4 4 17.0
2 38 4 18.6
3 72 4 24.0
4 106 4 20.0
5 140 4 20.0
6 174 4 15.6
7 208 4 16.0
8 242 4 13.0
9 276 4 8.2
10 310 4 8.0
1 4 5 23.0
2 38 5 21.0
3 72 5 18.0
4 106 5 17.0
5 140 5 16.2
6 174 5 14.0
7 208 5 14.2
8 242 5 13.4
9 276 5 11.8
10 310 5 11.4
6 4 6 10.4
7 38 6 12.3
8 72 6 13.2
9 106 6 11.6
10 140 6 8.4
4 4 7 22.8
5 38 7 22.4
6 72 7 21.4
7 106 7 18.8
8 140 7 18.3
9 174 7 16.2
10 208 7 15.0
1 4 8 22.2
2 38 8 20.0
3 72 8 21.0
4 106 8 23.0
5 140 8 16.8
6 174 8 11.0
7 208 8 13.0
8 242 8 17.0
9 276 8 13.0
10 310 8 12.6
1st col : class variable - hert test day(HTD)2nd col : fixed regression - days in milk(DIM)3rd col : animal4th col : test dayfat yield(TDY)

 

위 자료를 fixed_regression.dat로 저장

 

Pedigree

 

1 0 0
2 0 0
3 0 0
412
532
615
734
817
1st col : animal2nd col : sire3rd col : dam

 

위 혈통을 fixed_regression.ped로 저장

 

Left hand side와 Right hand side를 구하는 프로그램

 

PROGRAM fixed_regression_setup ! program name
! : fixed regression model setup program
! prgrammer
! : Park Byoungho
! date
! : 2010.1.27.

IMPLICIT NONE

! data dictionary INTEGER, PARAMETER :: no_of_cv = 1 ! number of class variables
INTEGER, PARAMETER :: no_of_rc = 5 ! number of regression coefficient
INTEGER, PARAMETER :: no_of_effect = no_of_cv + no_of_rc + 2 ! number of effects(2 = animal effect, permanent environment effect)

CHARACTER(LEN=256), PARAMETER :: data_filename = 'fixed_regression.dat' ! data file name
CHARACTER(LEN=256), PARAMETER :: pedi_filename = 'fixed_regression.ped' ! pedigree file name ! Legendre polynomial
REAL(KIND = 8), PARAMETER :: lambda(5,5) = &
RESHAPE((/SQRT(.5) , 0.0 , 0.0 , 0.0 , 0.0, &
0.0 , SQRT(1.5) , 0.0 , 0.0 , 0.0, &
-SQRT(2.5)*.5 , 0.0 , SQRT(2.5)*1.5 , 0.0 , 0.0, &
0.0 , -SQRT(3.5)*1.5, 0.0 , SQRT(3.5)*2.5, 0.0, &
SQRT(4.5)*3./8., 0.0 , -SQRT(4.5)*30./8., 0.0 , SQRT(4.5)*35./8./),(/5,5/)) !DO i = 1, 5
! WRITE(*,*) lambda(i,:)
!END DO

REAL(KIND = 8), ALLOCATABLE :: input_ind(:) ! storage for independent variable when reading data
REAL(KIND = 8), ALLOCATABLE :: ind(:) ! storage for independent variable after making Legendre polynimal
REAL(KIND = 8) :: dep ! storage for dependent variable when reading data INTEGER :: levels_of_cv ! levels of calss variable
REAL(KIND = 8) :: dim_max, dim_min, dim_range ! maximum, minimum and range of days in milk
REAL(KIND = 8) :: sut ! standardized unit of time
REAL(KIND = 8) :: sut_poly(5) ! polynomial of standardized unit of time
INTEGER :: no_of_animal ! number of animals
INTEGER :: no_of_eq ! number of equations REAL(KIND = 8), ALLOCATABLE :: rhs(:) ! right-hand side
INTEGER, ALLOCATABLE :: eq_no(:) ! storage for equation number REAL(KIND = 8) :: temp_ind ! temporary independent variable
INTEGER :: temp_eq_no INTEGER :: animal, sire, dam ! animal, sire, dam name for pedigree file INTEGER :: i,j ! for loop
INTEGER :: status ! I/O status
CHARACTER(LEN = 40) :: error_msg ! error message

ALLOCATE(input_ind(no_of_cv + 2)) ! independent variable
ALLOCATE(ind( no_of_cv + no_of_rc + 2)) ! independent variable after making Legendre polynimal
ALLOCATE(eq_no(no_of_cv + no_of_rc + 2)) ! storage for equation number ! open data file
OPEN(UNIT = 10, FILE = data_filename, STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg)

IF (status /= 0) THEN ! file open failed
WRITE (*,'(1X, A, A)') 'Data file open failed -- error message : ', error_msg
STOP
END IF ! look for the maximum and minimum of days in milk
! and levels of class variable(herd-test day)
! read first line of data file
READ(UNIT = 10, FMT = *, IOSTAT = status) input_ind, dep

levels_of_cv = input_ind(1)
dim_max = input_ind(2)
dim_min = input_ind(2)
no_of_animal = input_ind(3)

REWIND(10) ! read each line of data file
DO
READ(UNIT = 10, FMT = *, IOSTAT = status) input_ind
IF (status /= 0) EXIT ! end of file IF (input_ind(1) > levels_of_cv) THEN
levels_of_cv = input_ind(1)
END IF IF (input_ind(2) < dim_min) THEN
dim_min = input_ind(2)
END IF IF (input_ind(2) > dim_max) THEN
dim_max = input_ind(2)
END IF IF (input_ind(3) > no_of_animal) THEN
no_of_animal = input_ind(3)
END IF END DO REWIND(10) no_of_eq = levels_of_cv + no_of_rc + no_of_animal * 2
ALLOCATE(rhs(no_of_eq)) ! right-hand side

! initialize
rhs = 0.0 ! calculate the range of days in milk
dim_range = dim_max - dim_min

WRITE(*,*) 'Levels of class variable = ', levels_of_cv
WRITE(*,*) 'Minimum value of days in milk = ', dim_min
WRITE(*,*) 'Maximum value of days in milk = ', dim_max
WRITE(*,*) 'Range of days in milk = ', dim_range
WRITE(*,*) 'Number of animals = ', no_of_animal ! open file for fixed_regression_setup result
OPEN(UNIT = 15, FILE = 'fixed_regression_setup.rst', STATUS = 'REPLACE', ACTION = 'WRITE') ! write to fixed_regression_setup.rst
WRITE(15,*) 'Levels of class variable = ', levels_of_cv
WRITE(15,*) 'Minimum value of days in milk = ', dim_min
WRITE(15,*) 'Maximum value of days in milk = ', dim_max
WRITE(15,*) 'Range of days in milk = ', dim_range
WRITE(15,*) 'Number of animals = ', no_of_animal CLOSE(15) ! close

! left hand side
! open data file for writing left-hand side
OPEN(UNIT = 20, FILE = 'lhs.dat', STATUS = 'REPLACE', ACTION = 'WRITE')
! read each line
DO
READ(UNIT = 10, FMT = *, IOSTAT = status) input_ind, dep
IF (status /= 0) EXIT ! end of file ! assign equation number : class variable
eq_no(1) = input_ind(1) ! assign equation number : regression coefficient
DO i = 2, no_of_cv + no_of_rc
eq_no(i) = levels_of_cv + i - 1
END DO

! assign equation number : animal effect, permanent environmental effect
eq_no(no_of_cv + no_of_rc + 1) = levels_of_cv + no_of_rc + input_ind(3)
eq_no(no_of_cv + no_of_rc + 2) = levels_of_cv + no_of_rc + no_of_animal + input_ind(3) ! assign value : class variable
ind(1) = 1 ! assign value : regression coefficients
sut = -1. + 2. * (input_ind(2) - dim_min) / dim_range ! standardized unit of time
sut_poly(1) = 1.0 ! polynomial of standardized unit of time
sut_poly(2) = sut
sut_poly(3) = sut_poly(2) * sut
sut_poly(4) = sut_poly(3) * sut
sut_poly(5) = sut_poly(4) * sut
ind(2:6) = MATMUL(sut_poly, lambda)

! assign value : animal effect, permanent environmental effect
ind(7:8) = 1.0

! for fixed_effect ( X'X, Q'Q, Z'Z)
DO i = 1, no_of_cv + no_of_rc + 2
! swap equation number
temp_eq_no = eq_no(i)
eq_no(i) = eq_no(1)
eq_no(1) = temp_eq_no

! swap independent variable
temp_ind = ind(i)
ind(i) = ind(1)
ind(1) = temp_ind

WRITE(20,*) 1, (eq_no(j),ind(1)*ind(j), j=1, no_of_cv + no_of_rc + 2)

rhs(eq_no(1)) = rhs(eq_no(1)) + ind(1) * dep
END DO END DO ! for processing variance ratio of permanent environmental effect
DO i = levels_of_cv + no_of_rc + no_of_animal + 1, no_of_eq
WRITE(20,*) 2, i, (0, j = 3, 1 + no_of_effect * 2)
END DO CLOSE(10) ! close data file

! open pedigree file
OPEN(UNIT = 30, FILE = pedi_filename, STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg)

IF (status /= 0) THEN ! file open failed
WRITE (*,'(1X, A, A)') 'pedigree file open failed -- error message : ', error_msg
STOP
END IF ! read each line
DO
READ(UNIT = 30, FMT = *, IOSTAT = status) animal, sire, dam
IF (status /= 0) EXIT ! end of file ! find location of equations
animal = animal + levels_of_cv + no_of_rc
IF (sire /= 0) sire = sire + levels_of_cv + no_of_rc ! if sire is 0, sire is unknown
IF (dam /= 0) dam = dam + levels_of_cv + no_of_rc ! if dam is 0, dam is unknown ! write animal effect
IF ( sire /= 0 .AND. dam /= 0) THEN ! both parents are known
WRITE(UNIT = 20, FMT = *) 1001, animal, sire, dam, (0, i = 4, 1 + no_of_effect * 2)
WRITE(UNIT = 20, FMT = *) 1002, sire, animal, dam, (0, i = 4, 1 + no_of_effect * 2)
WRITE(UNIT = 20, FMT = *) 1002, dam, animal, sire, (0, i = 4, 1 + no_of_effect * 2)
ELSE IF (sire /= 0 .AND. dam == 0) THEN ! only sire is known
WRITE(UNIT = 20, FMT = *) 1003, animal, sire, dam, (0, i = 4, 1 + no_of_effect * 2)
WRITE(UNIT = 20, FMT = *) 1004, sire, animal, dam, (0, i = 4, 1 + no_of_effect * 2)
ELSE IF (sire == 0 .AND. dam /= 0) THEN ! only dam is known
WRITE(UNIT = 20, FMT = *) 1003, animal, dam, sire, (0, i = 4, 1 + no_of_effect * 2)
WRITE(UNIT = 20, FMT = *) 1004, dam, animal, sire, (0, i = 4, 1 + no_of_effect * 2)
ELSE ! no parents are known
WRITE(UNIT = 20, FMT = *) 1005, animal, sire, dam, (0, i = 4, 1 + no_of_effect * 2)
END IF
END DO ! end of reading data CLOSE(30) ! close animal pedigree file CLOSE(20) ! close left hand side

! open file for right hand side
OPEN(UNIT = 40, FILE = 'rhs.dat', STATUS = 'REPLACE', ACTION = 'WRITE') ! write right-hand side
DO i = 1, no_of_eq
WRITE(40, *) rhs(i)
END DO CLOSE(40) ! close END PROGRAM fixed_regression_setup

 

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

 

컴파일 : gfortran fixed_regression_setup.f95 -o fixed_regression_setup.exe
fixed_regression_setup.f95로 구한LHS 및 RHS

 

lhs.dat 참조 rhs.dat 참조

 

lhs.dat 의 정렬(유닉스 sort 명령어 사용)

 

sort -nk 2 lhs.dat -o sorted_lhs.dat

 

해를 구하는 프로그램(solver)

 

PROGRAM fixed_regression_solve ! program name
! : iteration on data of fixed regression model
! : solve program
! prgrammer
! : Park Byoungho
! date
! : 2010.1.27.

IMPLICIT NONE

! data dictionary INTEGER, PARAMETER :: no_of_cv = 1 ! number of class variables
INTEGER, PARAMETER :: no_of_rc = 5 ! number of regression coefficient
INTEGER, PARAMETER :: no_of_effect = no_of_cv + no_of_rc + 2 ! number of effects

REAL(KIND = 8), PARAMETER :: animal_vr = 3.710 / 5.521 ! variance ratio of animal genetic variance and residual variance
REAL(KIND = 8), PARAMETER :: pe_vr = 3.710 / 8.470 ! variance ratio of permanent environmental efeect variance and residual variance INTEGER, PARAMETER :: levels_of_cv = 10 ! levels of each calss variables
INTEGER, PARAMETER :: no_of_animal = 8
INTEGER, PARAMETER :: no_of_eq = levels_of_cv + no_of_rc + no_of_animal * 2 INTEGER :: no_of_lhs ! number of left hand side lines
REAL(KIND = 8), ALLOCATABLE :: lhs(:,:), rhs(:) ! left-hand side, right-hand side, x'y
REAL(KIND = 8), ALLOCATABLE :: solutions(:) ! solutions
INTEGER :: data_type ! fixed, animal, permanent environment INTEGER :: pre_eq_no ! previous equation number
INTEGER :: cur_eq_no ! current equation number REAL(KIND = 8) :: diag_ele ! diagonal element
REAL(KIND = 8) :: temp_rhs ! temporary rihgt hand side
REAL(KIND = 8) :: pre_sol ! previous solution INTEGER :: iteration ! iteration
REAL(KIND = 8) :: epsilon ! sum of squares of difference between old and new solutions
INTEGER :: i,j ! loop
INTEGER :: status ! i/o status
CHARACTER(LEN = 40) :: error_msg ! error message
INTEGER, PARAMETER :: MAX_ITER = 1000 ! maximum number of iteration
REAL(KIND = 8), PARAMETER :: criteria = 1.E-12 ! criteria for stopping ALLOCATE(rhs(no_of_eq))
ALLOCATE(solutions(no_of_eq))

solutions = 0 !open left-hand side file
OPEN(UNIT = 10, FILE = 'sorted_lhs.dat', STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg)
IF (status /= 0) THEN ! file open failed
WRITE (*,'(1X, A, A)') 'Sorted lhs file open failed -- error message : ', error_msg
STOP
END IF ! count the lines of lhs
no_of_lhs = 0
DO
READ(UNIT = 10, FMT = *, IOSTAT = status)
IF (status /= 0) EXIT ! reach end of file
no_of_lhs = no_of_lhs + 1
END DO ! decide array size of lhs
ALLOCATE(lhs(no_of_lhs, 1 + (no_of_cv + no_of_rc + 2) * 2))

!write(*,*) 'number of lines of left hand side = ', no_of_lhs ! store lhs to array
REWIND(10)
DO i = 1, no_of_lhs
READ(UNIT = 10, FMT = *, IOSTAT = status) lhs(i,:)
IF (status /= 0) EXIT ! reach end of file
END DO
CLOSE(10) !DO i = 1, no_of_lhs
! write(*,*) i, lhs(i,:)
!END DO !open right-hand side file
OPEN(UNIT = 20, FILE = 'rhs.dat', STATUS = 'OLD', ACTION = 'READ', IOSTAT = status, IOMSG = error_msg) IF (status /= 0) THEN ! file open failed
WRITE (*,'(1X, A, A)') 'RHS file open failed -- error message : ', error_msg
STOP
END IF

! read and store right-hand side
DO i = 1, no_of_eq
READ(UNIT = 20,FMT = *) rhs(i)
END DO
close(20)! write(*,*) 'right hand side ...'
! DO i = 1, levels_of_cv + no_of_rc
! write(*,*) i, rhs(i)
! END DO DO iteration = 1, MAX_ITER

pre_eq_no = 1
diag_ele = 0.0
temp_rhs = rhs(1)
epsilon = 0.0 ! start reading and processing each line of lhs
DO i = 1, no_of_lhs data_type = INT(lhs(i,1))
cur_eq_no = INT(lhs(i,2)) IF(pre_eq_no /= cur_eq_no) THEN
pre_sol = solutions(pre_eq_no) ! store old solution
solutions(pre_eq_no) = temp_rhs / diag_ele ! 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 diag_ele = 0.0
temp_rhs = rhs(cur_eq_no)
END IF ! process off-diagnoal
SELECT CASE(data_type)
CASE(1) ! fixed effect
! adjust right-hand side
! temp_rhs - independend variable * solutions
DO j = 2, no_of_effect
temp_rhs = temp_rhs - lhs(i, 2*j+1) * solutions(INT(lhs(i,2*j)))
END DO ! add diagonal element
diag_ele = diag_ele + lhs(i,3)

pre_eq_no = cur_eq_no
CASE(2) ! permanent environmental effect
diag_ele = diag_ele + pe_vr
pre_eq_no = cur_eq_no
CASE(1001) ! animal effect, both parents are known, animal is diagonal
temp_rhs = temp_rhs - solutions(INT(lhs(i,3))) * (-0.5) * 2.0 * animal_vr
temp_rhs = temp_rhs - solutions(INT(lhs(i,4))) * (-0.5) * 2.0 * animal_vr
diag_ele = diag_ele + 2.0 * animal_vr
pre_eq_no = cur_eq_no
CASE(1002) ! animal effect, both parents are known, sire or dam is diagonal
temp_rhs = temp_rhs - solutions(INT(lhs(i,3))) * (-0.5) * 2.0 * animal_vr
temp_rhs = temp_rhs - solutions(INT(lhs(i,4))) * 0.25 * 2.0 * animal_vr
diag_ele = diag_ele + 0.25 * 2.0 * animal_vr
pre_eq_no = cur_eq_no
CASE(1003) ! animal effect, one parent is known, animal is diagonal
temp_rhs = temp_rhs - solutions(INT(lhs(i,3))) * (-0.5) * (4.0/3.0) * animal_vr
diag_ele = diag_ele + (4.0/3.0) * animal_vr
pre_eq_no = cur_eq_no
CASE(1004) ! animal effect, one parents is known, parent is diagonal
temp_rhs = temp_rhs - solutions(INT(lhs(i,3))) * (-0.5) * (4.0/3.0) * animal_vr
diag_ele = diag_ele + 0.25 * (4.0/3.0) * animal_vr
pre_eq_no = cur_eq_no
CASE(1005) ! animal effect, no parents is known
diag_ele = diag_ele + animal_vr
pre_eq_no = cur_eq_no
END SELECT END DO
! end reading lhs ! calculate last solution
pre_sol = solutions(pre_eq_no) ! store old solution
solutions(pre_eq_no) = temp_rhs / diag_ele ! 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(*,*) iteration,'th iteration''s solutions'
! DO i = 1, no_of_rc
! WRITE(*,*) i, solutions(i)
! END DO epsilon = epsilon / (no_of_eq) ! write iteration number and epsilon
write(*,*) 'iteration = ', iteration , ', epsilon = ', epsilon
IF (epsilon < criteria) THEN
EXIT
END IF END DO
! end iteration ! open file for writing solution
OPEN(UNIT=30, FILE='sol.dat', STATUS='REPLACE', ACTION='WRITE', IOSTAT=status)
DO i = 1, no_of_eq
WRITE(30,*) i, solutions(i)
END DO
CLOSE(30)END PROGRAM fixed_regression_solve

 

위 프로그램을 fixed_regression_solve로 저장

 

컴파일 : gfortran fixed_regression_solve.f95 -o fixed_regression_solve.exe
프로그램 컴파일 및 실행화면

 

 

 

파일

1264520450_fixed_regression_model.zip

다운로드

1264520450_fixed_regression_model.zip