Fortran 95 source code to calculate a Kendall's Rank Correlation Coefficient.

calc_rcc_kendall.f95

!****************************************************
! ケンドールの順位相関係数の計算
!
! DATE          AUTHOR          VERSION
! 2019.12.13    mk-mode.com     1.00 新規作成
!
! Copyright(C) 2019 mk-mode.com All Rights Reserved.
!****************************************************
!
module cst
  ! SP: 単精度(4), DP: 倍精度(8)
  integer,     parameter :: SP = kind(1.0)
  integer(SP), parameter :: DP = selected_real_kind(2 * precision(1.0_SP))
end module cst

module comp
  use cst
  implicit none
  private
  public :: rcc_kendall

contains
  ! ケンドールの順位相関係数の計算
  !
  ! :param(in) real(8)  xs
  ! :param(in) real(8)  ys
  ! :return    real(8) rcc
  function rcc_kendall(xs, ys) result(rcc)
    use cst
    implicit none
    real(DP), intent(in) :: xs(:), ys(:)
    real(DP)    :: rcc
    integer(SP) :: n, p, q
    real(DP)    :: s_x, s_y, n2
    integer(SP), allocatable :: rank_x(:), rank_y(:)
    integer(SP), allocatable :: tai_x(:, :), tai_y(:, :)

    rcc = -1.0_DP
    n = size(xs)
    if (size(ys) /= n) return
    allocate(rank_x(n))
    allocate(rank_y(n))
    allocate(tai_x(n, 2))
    allocate(tai_y(n, 2))
    ! ランク付け
    call rank(xs, rank_x)
    call rank(ys, rank_y)
    ! P（x_s と x_t, y_s と y_t の大小関係が一致する組の数）
    ! Q（x_s と x_t, y_s と y_t の大小関係が不一致の組の数）
    ! （x_s = x_t or y_s = y_t は除く）
    call calc_pq(rank_x, rank_y, p, q)
    ! 同順位の数（インデックスが順位）
    call tai(rank_x, tai_x)
    call tai(rank_y, tai_y)
    ! Tx, Ty の sum 部分
    call sum_t(tai_x, s_x)
    call sum_t(tai_y, s_y)
    ! 相関係数
    n2 = (n * n - n) / 2.0_DP
    rcc = (p - q) / (sqrt(n2 - s_x) * sqrt(n2 - s_y))

    deallocate(rank_x)
    deallocate(rank_y)
    deallocate(tai_x)
    deallocate(tai_y)
  end function rcc_kendall

  ! ランク付け
  ! （同順位を中央（平均）順位(mid-rank)にする必要はない）
  !
  ! :param(in)  real(8)    src(:)
  ! :param(out) integer(4) res(:)
  subroutine rank(src, res)
    implicit none
    real(DP),    intent(in)  :: src(:)
    integer(SP), intent(out) :: res(:)
    integer(SP) :: c, i, j, s

    res = 0.0_DP
    s = size(src)
    do i = 1, s
      c = 0
      do j = 1, s
        if (src(i) < src(j)) c = c + 1
      end do
      res(i) = c + 1
    end do
  end subroutine rank

  ! P（x_s と x_t, y_s と y_t の大小関係が一致する組の数）
  ! Q（x_s と x_t, y_s と y_t の大小関係が不一致の組の数）
  ! （x_s = x_t or y_s = y_t は除く）
  subroutine calc_pq(rank_x, rank_y, p, q)
    implicit none
    integer(SP), intent(in)  :: rank_x(:), rank_y(:)
    integer(SP), intent(out) :: p, q
    integer(SP) :: i, j, n, w

    n = size(rank_x)
    p = 0
    q = 0
    do i = 1, n - 1
      do j = i + 1, n
        w = (rank_x(i) - rank_x(j)) * (rank_y(i) - rank_y(j))
        if (w > 0) then
          p = p + 1
        else if (w < 0) then
          q = q + 1
        end if
      end do
    end do
  end subroutine calc_pq

  ! 同順位の数
  !
  ! :param(in)  real(8)       src(:)
  ! :param(out) integer(r) res(:, 2)
  subroutine tai(src, res)
    implicit none
    integer(SP), intent(in)  :: src(:)
    integer(SP), intent(out) :: res(:, :)
    integer(SP) :: c, i, j, s
    integer(SP), allocatable :: wk(:)

    s = size(src)
    allocate(wk(s))
    wk = src
    res = 0
    do i = 1, s
      c = 0
      do j = i, s
        if (wk(j) == 0) cycle
        if (wk(i) == wk(j)) then
          c = c + 1
          if (i /= j) wk(j) = 0
        end if
      end do
      res(i, :) = (/wk(i), c/)
    end do
    deallocate(wk)
  end subroutine tai

  ! Tx, Ty の sum 部分
  !
  ! :param(in)  integer(4) src(:, 2)
  ! :param(out) real(8)          s_t
  subroutine sum_t(src, s_t)
    implicit none
    integer(SP), intent(in)  :: src(:, :)
    real(DP),    intent(out) :: s_t
    integer(SP) :: i, s
    integer(DP) :: t

    s = size(src) / 2
    s_t = 0.0_DP
    do i = 1, s
      t = src(i, 2)
      s_t = s_t + (t * t * t - t) / 2.0_DP
    end do
  endsubroutine sum_t
end module comp

program calc_rcc_kendall
  use cst
  use comp
  implicit none
  real(DP)    :: rcc
  ! === タイ（同順位）が存在しない例
  !real(DP), parameter :: X(10) = (/ &
  !  & 1.0_DP, 2.0_DP, 3.0_DP, 4.0_DP,  5.0_DP, &
  !  & 6.0_DP, 7.0_DP, 8.0_DP, 9.0_DP, 10.0_DP  &
  !& /)
  !real(DP), parameter :: Y(10) = (/ &
  !  & 1.0_DP, 3.0_DP, 5.0_DP, 7.0_DP,  9.0_DP, &
  !  & 2.0_DP, 4.0_DP, 6.0_DP, 8.0_DP, 10.0_DP  &
  !& /)

  ! === タイ（同順位）が存在する例
  real(DP), parameter :: X(10) = (/ &
    & 1.0_DP, 2.0_DP, 3.0_DP, 4.0_DP,  5.0_DP, &
    & 5.0_DP, 7.0_DP, 8.0_DP, 9.0_DP, 10.0_DP  &
  & /)
  real(DP), parameter :: Y(10) = (/ &
    & 1.0_DP, 3.0_DP, 5.0_DP, 6.0_DP,  9.0_DP, &
    & 2.0_DP, 4.0_DP, 6.0_DP, 8.0_DP, 10.0_DP  &
  & /)

  ! === サイズが異なる例 ( => 実行時にエラー )
  !real(DP), parameter :: X(10) = (/ &
  !  & 1.0_DP, 2.0_DP, 3.0_DP, 4.0_DP,  5.0_DP, &
  !  & 6.0_DP, 7.0_DP, 8.0_DP, 9.0_DP, 10.0_DP  &
  !& /)
  !real(DP), parameter :: Y(9) = (/ &
  !  & 1.0_DP, 3.0_DP, 5.0_DP, 7.0_DP,  9.0_DP, &
  !  & 2.0_DP, 4.0_DP, 6.0_DP, 8.0_DP           &
  !& /)

  ! === X のサイズがゼロの例（ => 当然、コンパイルエラー ）
  !real(DP), parameter :: X(0) = (//)
  !real(DP), parameter :: Y(9) = (/ &
  !  & 1.0_DP, 3.0_DP, 5.0_DP, 7.0_DP,  9.0_DP, &
  !  & 2.0_DP, 4.0_DP, 6.0_DP, 8.0_DP, 10.0_DP  &
  !& /)

  ! === 数値以外のものが存在する例（ => 当然、コンパイルエラー ）
  !real(DP), parameter :: X(10) = (/ &
  !  & 1.0_DP, 2.0_DP, 3.0_DP, 4.0_DP,  5.0_DP, &
  !  & 6.0_DP, 7.0_DP, 8.0_DP, 9.0_DP, 10.0_DP  &
  !& /)
  !real(DP), parameter :: Y(10) = (/ &
  !  & 1.0_DP, 3.0_DP, 5.0_DP, 7.0_DP,  9.0_DP, &
  !  & "ABC", 4.0_DP, 6.0_DP, 8.0_DP, 10.0_DP  &
  !& /)

  rcc = rcc_kendall(X, Y)
  print '("X = ", 10F5.1)', X
  print '("Y = ", 10F5.1)', Y
  if (rcc == -1.0_DP) then
    print *, "[ERROR]"
  else
    print '("Spearman''s RCC = ", F11.8)', rcc
  end if
end program calc_rcc_kendall