日本畜産統計戦略会議（JCSAS）

変量効果の推定とBLUP法 単行本 2007/2
佐々木 義之 (著)
ISBN-13: 978-4876987023

のp124にて掲載された内容を実装したFortran90プログラム部品です。自由にご利用ください。

---

!tm -------------------------------------------------------
  subroutine create_a00
! numerator relatonship matrix subroutine by Fortran90
! 出典：
! 変量効果の推定とBLUP法 単行本 2007/2
! 佐々木 義之 (著)
! ISBN-13: 978-4876987023
! p124
!----------------------------
! a s d sex   -> a s d renumbered!
! 1 0 0 m        1 0 0
! 2 1 0 m        2 1 0
! x 1 0 f        3 1 0
! y 1 0 f        4 1 0
! 3 1 0 m        5 2 0
! 4 1 0 m        6 1 0
! 5 3 0 m        7 5 0
! 6 0 x m        8 0 3
! 7 3 y m        9 5 4
!----------------------------
    integer, parameter :: nanim = 9
    integer :: i, j
    real (8) :: amat(0:nanim, 0:nanim)
    real (8) :: atmp
    integer :: sid(nanim), did(nanim)
    print *, '--- create_A00_sasaki ---'
!tm -------------------------
    sid(1) = 0
    sid(2) = 1
    sid(3) = 1
    sid(4) = 1
    sid(5) = 2
    sid(6) = 1
    sid(7) = 5
    sid(8) = 0
    sid(9) = 5
!tm -------------------------
    did(1) = 0
    did(2) = 0
    did(3) = 0
    did(4) = 0
    did(5) = 0
    did(6) = 0
    did(7) = 0
    did(8) = 3
    did(9) = 4
!tm -------------------------
    do i = 1, nanim
      if (sid(i)>0 .and. did(i)>0) then
        amat(i, i) = 1.0d0 + 0.5d0*amat(sid(i), did(i))
      else
        amat(i, i) = 1.0d0
      end if
      do j = 1, i - 1
        atmp = 0.0d0
        if (sid(i)>0) atmp = atmp + 0.5d0*amat(j, sid(i))
        if (did(i)>0) atmp = atmp + 0.5d0*amat(j, did(i))
        amat(j, i) = atmp
        amat(i, j) = atmp
      end do
    end do
    write (*, *) '2016 blup sasaki book p124'
    write (*, *) '* a s d sex -> a s d revnum'
    write (*, *) '* 1 0 0 m      1 0 0'
    write (*, *) '* 2 1 0 m      2 1 0'
    write (*, *) '* x 1 0 f      3 1 0'
    write (*, *) '* y 1 0 f      4 1 0'
    write (*, *) '* 3 1 0 m      5 2 0'
    write (*, *) '* 4 1 0 m      6 1 0'
    write (*, *) '* 5 3 0 m      7 5 0'
    write (*, *) '* 6 0 x m      8 0 3'
    write (*, *) '* 7 3 y m      9 5 4'
    write (*, *) ''
    do i = 1, nanim
      write (*, '(20f7.4)')(amat(i,j), j=1, nanim)
    end do
! 1.0000 0.5000 0.5000 0.5000 0.2500 0.5000 0.1250 0.2500 0.3750
! 0.5000 1.0000 0.2500 0.2500 0.5000 0.2500 0.2500 0.1250 0.3750
! 0.5000 0.2500 1.0000 0.2500 0.1250 0.2500 0.0625 0.5000 0.1875
! 0.5000 0.2500 0.2500 1.0000 0.1250 0.2500 0.0625 0.1250 0.5625
! 0.2500 0.5000 0.1250 0.1250 1.0000 0.1250 0.5000 0.0625 0.5625
! 0.5000 0.2500 0.2500 0.2500 0.1250 1.0000 0.0625 0.1250 0.1875
! 0.1250 0.2500 0.0625 0.0625 0.5000 0.0625 1.0000 0.0313 0.2813
! 0.2500 0.1250 0.5000 0.1250 0.0625 0.1250 0.0313 1.0000 0.0938
! 0.3750 0.3750 0.1875 0.5625 0.5625 0.1875 0.2813 0.0938 1.0625
  end subroutine create_a00

---