Fortran - ケンドール順位相関係数の計算!

Updated:


Fortran 95 でケンドールの順位相関係数(Kendall’s Rank Correlation Coefficient)の計算をしてみました。

0. 前提条件

  • Debian GNU/Linux 10.3 (64bit) での作業を想定。
  • GCC 9.2.0 (GFortran 9.2.0) でのコンパイルを想定。

1. ケンドールの順位相関係数について

\(n\) 対の順位データ \((x_i, y_i) (i=1,2,\cdots,n)\) があるとき、それの中から取り出した \((x_s, y_s),\ (x_t, y_t)\ (1 \leq s \lt t \leq n)\) において、

\[\begin{eqnarray*} P&:& x_s と x_t,\ y_s と y_t の大小関係が一致する組の数 \\ Q&:& x_s と x_t,\ y_s と y_t の大小関係が不一致の組の数 \\ N&:& 組の総数=\binom{n}{2}=\displaystyle _nC_2=\frac{n(n-1)}{2} \end{eqnarray*}\]

とおくと、 ケンドールの順位相関係数(Kendall’s Rank Correlation Coefficient) \(\tau_a, \tau_b\) は次式で表される。(ケンドールの \(\tau_a\) (Kendall’s tau a), ケンドールの \(\tau_b\) (Kendall’s tau b)とも呼ばれる)

(1)同順位(タイ)が存在しない場合、 \(\begin{eqnarray*} \tau_a = \frac{P - Q}{N} \end{eqnarray*}\)

(2)同順位(タイ)が存在する場合、 \(\begin{eqnarray*} \tau_b &=& \frac{P - Q}{\sqrt{N - T_x}\sqrt{N - T_y}} \\ 但し、\ \ T_x &=& \displaystyle \sum_{i=1}^{n_x}\frac{t_i({t_i}^2 - 1)}{2} \\ T_y &=& \displaystyle \sum_{j=1}^{n_y}\frac{t_j({t_j}^2 - 1)}{2} \\ &&(n_x,\ n_y\ は同順位の数、\ t_i,\ t_j\ は同順位となる順位の個数) \end{eqnarray*}\)

また、次式で表されるものを グッドマン=クラスカルの \(\gamma\) (Goodman and Kruskal’s gamma) と呼ぶ。

\[\begin{eqnarray*} \gamma = \frac{P - Q}{P + Q} \end{eqnarray*}\]

(注1)同順位(タイ)がない場合、 \(T_x = T_y = 0,\ P + Q = N\) となり、結果として、 \(\tau_a = \tau_b = \gamma\) となる。
(注2)\(\tau_c\) なる式もあるが、 \(\tau_a,\ \tau_b,\ \gamma\) とやや性質が異なるので今回は割愛。

2. Fortran ソースコードの作成

  • 計算の本質は comp モジュール。
  • 使用する2つの配列 X, Y を適宜変更する等する。(calc_rcc_kendall プログラム内)

File: calc_rcc_kendall.f95

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
!****************************************************
! ケンドールの順位相関係数の計算
!
! 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

3. コンパイル&リンク

$ gfortran -Wall -O2 -o calc_rcc_kendall calc_rcc_kendall.f95

4. プログラムの実行

$ ./calc_rcc_kendall
X =   1.0  2.0  3.0  4.0  5.0  5.0  7.0  8.0  9.0 10.0
Y =   1.0  3.0  5.0  6.0  9.0  2.0  4.0  6.0  8.0 10.0
Spearman's RCC =  0.64285714

以上。





 

Sponsored Link

 

Comments