#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。
フィッシャーの正確確率検定
http://aoki2.si.gunma-u.ac.jp/R/fisher.html
ファイル名: fisher.jl 関数名: fisher
翻訳するときに書いたメモ
hybrid は省略した。
==========#
using Rmath, Printf
function fisher(x, y = NaN; exact = true, method = "Fisher") # or method = "Pearson"
function found()
if method == "Fisher"
nprod = temp - sum(perm_table[um .+ 1])
if nprod <= criterion + EPSILON
nntrue += exp(nprod - nntrue2 * log_expmax)
while nntrue >= EXPMAX
nntrue /= EXPMAX
nntrue2 += 1
end
end
else
hh = sum((um .- ex) .^ 2 ./ ex)
if hh >= stat_val || abs(hh - stat_val) <= EPSILON
nprod = temp - sum(perm_table[um .+ 1])
nntrue += exp(nprod - nntrue2 * log_expmax)
while nntrue >= EXPMAX
nntrue /= EXPMAX
nntrue2 += 1
end
end
end
ntables += 1
end
function search(x, y)
if y == 1
found()
elseif x == 1
t = um[1, 1] - um[y, 1]
if t >= 0
um[1, 1] = t
search(nc, y - 1)
um[1, 1] += um[y, 1]
end
else
search(x - 1, y)
while um[y, 1] != 0 && um[1, x] != 0
um[y, x] += 1
um[y, 1] -= 1
um[1, x] -= 1
search(x - 1, y)
end
um[y, 1] += um[y, x]
um[1, x] += um[y, x]
um[y, x] = 0
end
end
function exacttest()
denom2 = 0
denom = perm_table[n + 1] - sum(perm_table[ct .+ 1])
while denom > log_expmax
denom -= log_expmax
denom2 += 1
end
denom = exp(denom)
um[:, 1] = rt
um[1, :] = ct
search(nc, nr)
@printf("% s の方法による，正確な P 値 = % .10g\n", method, nntrue / denom * EXPMAX ^ (nntrue2 - denom2))
@printf("査察した分割表の個数は % s 個\n", ntables)
end
function chisqtest(t)
return sum((t .- ex) .^ 2 ./ ex)
end
function asymptotic()
@printf("カイ二乗値 = % g, 自由度 = % i, P 値 = % g\n", stat_val, df, pchisq(stat_val, df, false))
end
if ndims(x) == 2
t = x
else
indexy, indexx, t = table(y, x)
end
EPSILON = 1e-10
EXPMAX = 1e100
log_expmax = log(EXPMAX)
nr, nc = size(t)
um = zeros(Int, nr, nc)
rt = sum(t, dims=2)
ct = sum(t, dims=1)
n = sum(t)
ex = (rt * ct) ./ n
stat_val = chisqtest(t)
df = (nr - 1) * (nc - 1)
asymptotic()
if exact
perm_table = cumsum(vcat(0, log.(1:n + 1)))
temp = sum(perm_table[rt .+ 1])
criterion = temp - sum(perm_table[t .+ 1])
ntables = nntrue = nntrue2 = 0
exacttest()
end
end
function table(x, y) # 二次元
indicesx = sort(unique(x))
indicesy = sort(unique(y))
counts = zeros(Int, length(indicesx), length(indicesy))
for (i, j) in zip(indexin(x, indicesx), indexin(y, indicesy))
counts[i, j] += 1
end
return indicesx, indicesy, counts
end
分割表を与える場合
x = [5 3 2 1
4 3 5 2
2 3 1 2]
fisher(x)
# カイ二乗値 = 3.39631, 自由度 = 6, P 値 = 0.757711
# Fisher の方法による，正確な P 値 = 0.8091124268
# 査察した分割表の個数は 24871 個
fisher(x, method = "Pearson")
# カイ二乗値 = 3.39631, 自由度 = 6, P 値 = 0.757711
# Pearson の方法による，正確な P 値 = 0.7878188077
# 査察した分割表の個数は 24871 個
2 つのベクトルを与える場合
x = [1, 1, 1, 3, 2, 1, 1, 3, 3, 2, 3, 3, 1, 1, 3, 1, 1, 1, 2, 1,
1, 2, 2, 3, 1, 1, 2, 3, 3, 2];
y = [1, 2, 2, 3, 2, 1, 2, 2, 2, 3, 3, 2, 2, 1, 2, 1, 2, 2, 3, 1,
1, 1, 1, 3, 1, 2, 2, 3, 1, 3];
fisher(x, y)
# カイ二乗値 = 9.17504, 自由度 = 4, P 値 = 0.0568702
# Fisher の方法による，正確な P 値 = 0.03320873103
# 査察した分割表の個数は 1353 個
fisher(x, y, method="Pearson")
# カイ二乗値 = 9.17504, 自由度 = 4, P 値 = 0.0568702
# Pearson の方法による，正確な P 値 = 0.05323058543
# 査察した分割表の個数は 1353 個
