#==========
Julia の修行をするときに,いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。
一元配置分散分析(exact test)
http://aoki2.si.gunma-u.ac.jp/R/exact-oneway-test.html
ファイル名: exactonewaytest.jl 関数名: exactonewaytest
翻訳するときに書いたメモ
==========#
using Statistics, Rmath, Printf
function exactonewaytest(x, y=[]; permutation = true)
function found()
hh = performtest(um)
if hh <= pvalue + EPSILON
nprod = sum(perm_table[rt .+ 1]) - sum(perm_table[um .+ 1])
nntrue += exp(nprod - nntrue2 * log_expmax)
while nntrue >= EXPMAX
nntrue /= EXPMAX
nntrue2 += 1
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 permutationtest()
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)
pvalue = nntrue / denom * EXPMAX ^ (nntrue2 - denom2)
@printf("並べ替え検定による P 値 = % .10g\n", pvalue)
@printf("査察した分割表の個数は % s 個\n", ntables)
return pvalue
end
function performtest(u)
x = Array{Array{Float64,1},1}(undef, nr)
for i in 1:nr
x[i] = score[u[i, :] .> 0]
end
onewaytest(x)
end
function simpleonewaytest()
@printf("観察値による一元配置分散分析の P 値 = % g\n", pvalue)
end
if length(y) == 0
ni = map(length, x)
y = rep(1:length(ni), ni)
indexy, score, t = table(y, vcat(x...))
else
indexy, score, 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 = transpose(sum(t, dims=1))
n = sum(t)
q = cumsum(vcat(0, ct[1:nc-1])) .+ (ct .+ 1) .* 0.5
half = (n + 1) * 0.5
pvalue = performtest(t)
simpleonewaytest()
if permutation
perm_table = cumsum(vcat(0, log.(1:n + 1)))
ntables = nntrue = nntrue2 = 0
permutationtest()
end
end
function onewaytest(x)
n = map(length, x)
m = map(mean, x)
v = map(var, x)
w = n ./ v
sumw = sum(w)
tmp = sum((1 .- (w / sumw)) .^ 2 ./ (n .- 1)) / (nr ^ 2 - 1)
m0 = sum(w .* m) / sumw
F = sum(w .* (m .- m0) .^ 2) / ((nr - 1) * (1 + 2 * (nr - 2) * tmp))
df1 = nr - 1
df2 = 1 / (3 * tmp)
pf(F, df1, df2, false)
end
function rep(x, n::Array{Int64,1})
length(x) == length(n) || error("length(x) wasn't length(n)")
vcat([repeat([x[i]], n[i]) for i = 1:length(x)]...)
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 = [[36.7, 52.4, 65.8], [45.7, 61.9, 65.3], [52.6, 76.6, 81.3]];
exactonewaytest(x)
# 観察値による一元配置分散分析の P 値 = 0.440037
# 並べ替え検定による P 値 = 0.4107142857
# 査察した分割表の個数は 1680 個
※コメント投稿者のブログIDはブログ作成者のみに通知されます