#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。
Python による統計処理
http://aoki2.si.gunma-u.ac.jp/Python/PythonDir.tar.gz
一元配置分散分析（三群以上の平均値の差の検定）
http://aoki2.si.gunma-u.ac.jp/Python/oneway_ANOVA.pdf
ファイル名: onewayanova.jl
関数名: onewayanova(x, g; varequal=false)
onewayanova(ni, mi, ui, gi; varequal=false)
翻訳するときに書いたメモ
Python --> Julia は，R --> Julia より面倒。
==========#
using Statistics, Rmath, Printf
function onewayanova(x, g; varequal=false)
gi, ni = table(g)
k = length(ni)
z = Array{Array{Float64,1},1}(undef, k)
for (i, g2) in enumerate(gi)
z[i] = x[g .== g2]
end
mi = map(mean, z)
vi = map(var, z)
method = "One-way analysis of means"
if varequal
n = length(x)
mx = mean(x)
df1, df2 = k - 1, n - k
F = (sum(ni .* (mi .- mx).^2) / df1) / (sum((ni .- 1) .* vi) / df2)
else
method *= " (not assuming equal variances)"
wi = ni ./ vi
sum_wi = sum(wi)
tmp = sum((1 .- wi ./ sum_wi).^2 ./ (ni .- 1)) ./ (k^2 - 1)
m = sum(wi .* mi) / sum_wi
df1, df2 = k - 1, 1 / 3tmp
F = sum(wi .* (mi .- m).^2) / (df1 * (1 + 2 * (k - 2) * tmp))
end
output1(gi, ni, mi, vi)
output2(method, F, df1, df2)
end
function onewayanova(ni, mi, ui, gi; varequal=false)
x = []
g = []
for i in 1:length(ni)
sdi = sqrt(ui[i])
y = rnorm(ni[i], mi[i], sdi)
y = ((y .- mean(y)) / std(y)) .* sdi .+ mi[i]
append!(x, y)
append!(g, repeat([gi[i]], ni[i]))
end
onewayanova(x, g; varequal)
end
function table(x) # indices が少ないとき
indices = sort(unique(x))
counts = zeros(Int, length(indices))
for i in indexin(x, indices)
counts[i] += 1
end
return indices, counts
end
function output1(gi, ni, mi,vi)
println("Sumary of the data.")
@printf("%10s: %5s %15s %15s %15s\n",
"group", "n", "mean", "variance", "s.d.")
for i = 1:length(gi)
@printf("%10s: %5d %15.7f %15.7f %15.7f\n",
gi[i], ni[i], mi[i], vi[i], sqrt(vi[i]))
end
end
function output2(method, F, df1, df2)
println(method)
p = pf(F, df1, df2, false)
if df2 == floor(df2)
@printf("F = %.7g, d.f.1 = %d, d.f.2 = %d, p value = %.7g\n",
F, df1, df2, p)
else
@printf("F = %.7g, d.f.1 = %d, d.f.2 = %.7g, p value = %.7g\n",
F, df1, df2, p)
end
end
x = [2, 1, 2, 3, 4, 3, 5, 4, 6, 7, 2, 4];
g = ["a", "a", "a", "a", "b", "c", "c", "b", "b", "c", "c", "c"];
onewayanova(x, g, varequal=true)
# Sumary of the data.
# group: n mean variance s.d.
# a: 4 2.0000000 0.6666667 0.8164966
# b: 3 4.6666667 1.3333333 1.1547005
# c: 5 4.2000000 3.7000000 1.9235384
# One-way analysis of means
# F = 3.57149, d.f.1 = 2, d.f.2 = 9, p value = 0.0721381
onewayanova(x, g)
# Sumary of the data.
# group: n mean variance s.d.
# a: 4 2.0000000 0.6666667 0.8164966
# b: 3 4.6666667 1.3333333 1.1547005
# c: 5 4.2000000 3.7000000 1.9235384
# One-way analysis of means (not assuming equal variances)
# F = 6.256838, d.f.1 = 2, d.f.2 = 5.083311, p value = 0.04258988
onewayanova([10, 10, 10], [1.2, 3.2, 2.5], [1.1, 1.2, 1.3], ["a", "b", "c"], varequal=true)
# Sumary of the data.
# group: n mean variance s.d.
# a: 10 1.2000000 1.1000000 1.0488088
# b: 10 3.2000000 1.2000000 1.0954451
# c: 10 2.5000000 1.3000000 1.1401754
# One-way analysis of means
# F = 8.583333, d.f.1 = 2, d.f.2 = 27, p value = 0.001302067
onewayanova([10, 10, 10], [1.2, 3.2, 2.5], [1.1, 1.2, 1.3], ["a", "b", "c"])
# Sumary of the data.
# group: n mean variance s.d.
# a: 10 1.2000000 1.1000000 1.0488088
# b: 10 3.2000000 1.2000000 1.0954451
# c: 10 2.5000000 1.3000000 1.1401754
# One-way analysis of means (not assuming equal variances)
# F = 8.688278, d.f.1 = 2, d.f.2 = 17.97905, p value = 0.002290091
