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

calc_rcc_spearman.f95

!****************************************************
! スピアマンの順位相関係数の計算
!
! DATE          AUTHOR          VERSION
! 2019.12.12    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_spearman

contains
  ! スピアマンの順位相関係数の計算
  !
  ! :param(in) real(8)  xs
  ! :param(in) real(8)  ys
  ! :return    real(8) rcc
  function rcc_spearman(xs, ys) result(rcc)
    use cst
    implicit none
    real(DP), intent(in) :: xs(:), ys(:)
    real(DP)    :: rcc
    integer(SP) :: i, n, n3
    real(DP)    :: s_x, s_y, t_x, t_y, s
    real(DP), allocatable :: rank_x(:), rank_y(:)
    real(DP), 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)
    ! 同順位の数（インデックスが順位）
    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)
    ! 相関係数
    n3 = n * n * n - n
    t_x = (n3 - s_x) / 12.0_DP
    t_y = (n3 - s_y) / 12.0_DP
    s = 0.0_DP
    do i = 1, n
      s = s + (xs(i) - ys(i)) *  (xs(i) - ys(i))
    end do
    rcc = (t_x + t_y - s) / (2.0_DP * sqrt(t_x * t_y))

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

  ! ランク付け
  !
  ! :param(in)  real(8) src(:)
  ! :param(out) real(8) res(:)
  subroutine rank(src, res)
    implicit none
    real(DP), intent(in)  :: src(:)
    real(DP), intent(out) :: res(:)
    integer(SP) :: c, i, j, s
    real(DP), allocatable :: wk(:)

    s = size(src)
    allocate(wk(s))
    wk = 0.0_DP
    do i = 1, s
      c = 0
      do j = 1, s
        if (src(i) < src(j)) c = c + 1
      end do
      wk(i) = real(c + 1, DP)
    end do
    ! 同順位を中央（平均）順位(mid-rank)に
    do i = 1, s
      c = 0
      do j = 1, s
        if (wk(i) == wk(j)) c = c + 1
      end do
      res(i) = sum((/(j, j=int(wk(i)),int(wk(i))+c-1)/)) / real(c, DP)
    end do
    deallocate(wk)
  end subroutine rank

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

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

  ! Tx, Ty の sum 部分
  !
  ! :param(in)  real(8) src(:, 2)
  ! :param(out) real(8)       s_t
  subroutine sum_t(src, s_t)
    implicit none
    real(DP), 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 = int(src(i, 2), DP)
      if (t < 2) cycle
      s_t = s_t + t * t * t - t
    end do
  endsubroutine sum_t
end module comp

program calc_rcc_spearman
  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_spearman(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_spearman