スピアマンの順位相関係数の分布     Last modified: May 15, 2002

 今回は,うまく行かなかったものについて。

 今までのはうまい方法を見つけることができたので,短時間で目的を達することができましたが,スピアマンの順位相関係数はそのような方法を現在まで見つけられません。

 素直にこつこつと計算すれば時間はかかってもできないことはないのですが,難関は全ての順列を発生させることでした。これは,奥村晴彦著「アルゴリズム辞典」に載っていました。しかし,いずれにしろ,このアルゴリズムでは階乗に比例します。awk では ie = 7 が我慢の限界でしょう。


# spearman.awk

BEGIN {
   ie = 7 # 大きくしすぎないように
   for (n = 3; n <= ie; n++) {
       nn = (n^3-n)/3
       for (i = 0; i <= nn; i++) {
           t[i] = 0
       }
       genperm4()
       all = factorial(n)
       cum1 = 0
       cum2 = 1
       printf "\n\n Spearman Rank Correlation Coefficient    n = %i\n\n" \
              "       d^2      r           Freq.      p" \
              "       Cum. p     Cum. p'\n", n
       for (i = 0; i <= nn; i += 2) {
           p2 = t[i]/all
           cum1 += p2
           printf "%10i  %8.5f  %10li  %9.6f  %9.6f  %9.6f\n",
                   i, 1-6*i/(n^3-n), t[i], p2, cum1, cum2
           cum2 -= p2
       }
   }
}

# 順列を全て発生させ,d^2 を計算しその値別にカウントする
function genperm4(     i, i2, sum, k, tt, c)
{
   count = 0
   for (i = 0; i < n; i++) {
       c[i+1] = p[i] = i + 1
   }
   k = 1
   while (k < n) {
       if (k%2 == 1) {
           i = c[k]
       }
       else {
           i = 0
       }
       tt = p[k]
       p[k] = p[i]
       p[i] = tt
       d2 = 0
       for (i2 = 0; i2 < n; i2++) {
           d2 += (i2+1-p[i2])^2
       }
       t[d2]++
       k = 1
       while (c[k] == 0) {
           c[k] = k
           k++
       }
       c[k]--
   }
}

# 整数値 n(n ≧ 0 )の階乗 n! を求める
# 使用例: x = factorial(10)

function factorial(n,     retv, i)
{
   retv = 1
   for (i = 1; i <= n; i++) {
       retv *= i
   }
   return retv
}

・ 直前のページへ戻る  ・ E-mail to Shigenobu AOKI