Fisher の正確確率検定の拡張     Last modified: May 15, 2002

2 × 2 分割表において独立性の検定のために統計的検定手法に,Fisher の正確確率検定 というのがあります。以下のスクリプト(プログラム)は,K×L分割表への拡張版です。

この方面は,Mehta らが最近手を付けていますが,もっと古くからあります。このアルゴリズムは彼らのとは違います。

このスクリプトでも,再帰関数を使用しています。任意のサイズの分割表に対応するためにはこのようにします。

他の検定の場合には,chi_sq 関数の部分を作り替えればよいのです。

データ入力のしかたはシステムに合うように適当に変えてください。


# exact.awk

BEGIN {
  EPSILON = 1e-10 # 統計量の大小関係の判定のときに使う
  nntrue = loop = ntot = 0

  printf "1: Fisher の方法による  2: Pearson の方法による ... " > "/dev/stderr"
  getline Fisher < "/dev/stdin"

  kg = get_num("行数 ... ")
  kc = get_num("列数 ... ")

  for (j = 0; j < kg; j++) {
    for (i = 0; i < kc; i++) {
      prompt = sprintf("%i行%i列目の観察数 ... ", j+1, i+1)
      u[j,i] = temp = get_num(prompt)
      n_1_k[j] += temp
      k_tie[i] += temp
      ntot += temp
    }
  }
  for (j = 0; j < kg; j++) {
    for (i = 0; i < kc; i++) {
      expectation[j,i] = n_1_k[j]*k_tie[i]/ntot
      um[j,i] = 0
    }
  }
  for (i = 0; i < kc; i++) {
    um[0,i] = k_tie[i]
  }
  for (i = 1; i < kg; i++) {
    um[i,0] = n_1_k[i]
  }
  table[0] = 1
  for (i = 1; i <= ntot; i++) {
    table[i] = table[i-1]*i
  }
  criterion = calc_prob(u)
  chi_sq_v = chi_sq(u)
  denom = table[ntot]
  for (i = 0; i < kc; i++) {
    denom /= table[k_tie[i]]
  }
  kc1 = kc-1
  kg1 = kg-1
  gen_tab(kc1, kg1)
  str = "Sum of the binomial products of all the tables in "
  printf "Exact p = %#.15g ( = Sw / S ) %s\n\n", nntrue/denom, Fisher == 1 ? "Fisher" : "Pearson"
  printf "%sSw = %#.15g\n", str, nntrue
  printf "%sS  = %#.15g\n", str, denom
  printf "Number of tables generated: %i\n", loop
}

function calc_prob(um,     nprod, i, j)
{
  nprod = 1
  for (j = 0; j < kg; j++) {
    nprod *= table[n_1_k[j]]
    domi = table[um[j,0]]
    for (i = 1; i < kc; i++) {
      domi *= table[um[j,i]]
    }
    nprod /= domi
  }
  return nprod
}

function found(i, j,     nprod)
{
  if (Fisher == 1) {
  	nprod = calc_prob(um)
  	if (nprod <= criterion+EPSILON) {
  	  nntrue += nprod
  	}
  }
  else {
    if (chi_sq(um) >= chi_sq_v-EPSILON) {
      nprod = calc_prob(um)
      nntrue += nprod
    }
  }
  loop++
}

function chi_sq(u,    i, j, temp, retv)
{
  retv = 0
  for (j = 0; j < kg; j++) {
    for (i = 0; i < kc; i++) {
      temp = expectation[j,i] - u[j,i]
      retv += temp*temp / expectation[j,i]
    }
  }
  return retv
}

function gen_tab(x, y,     t)
{
  if (y == 0) {
    found()
  }
  else if (x == 0) {
    t = um[0,0]-um[y,0]
    if (t < 0) {
      return
    }
    um[0,0] = t
    gen_tab(kc1, y-1)
    um[0,0] += um[y,0]
  }
  else {
    gen_tab(x-1, y)
    while (um[y,0] && um[0,x]) {
      um[y,x]++
      um[y,0]--
      um[0,x]--
      gen_tab(x-1, y)
    }
    um[y,0] += um[y,x]
    um[0,x] += um[y,x]
    um[y,x] = 0
  }
}

function get_num(prompt,   val)
{
  printf prompt > "/dev/stderr"
  getline val < "/dev/stdin"
  return val
}


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