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
}