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 }