#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。
中央値の（差の）信頼区間
http://aoki2.si.gunma-u.ac.jp/R/CI4median.html
ファイル名: ci4median.jl 関数名: ci4median
翻訳するときに書いたメモ
Rmath が使えることがわかったのは大収穫だった。
==========#
using Distributions, Rmath
function ci4median(x, y; conflevel = 0.95,
method = "independent.sample" #= or "paired.sample" =#)
ci4median(x; y=y, conflevel, method)
end
function ci4median(x; y=NaN, conflevel = 0.95, method = "independent.sample")
0 < conflevel < 1 || error("0 < conflevel < 1")
length(y) == 1 && (method = "one.sample")
method in ["independent.sample", "paired.sample", "one.sample"] || error("check 'method'")
if method == "independent.sample"
n1 = length(x)
n2 = length(y)
vector = sort!(vec(x .- y'));
k = n1 <= 100 && n2 <= 100 ?
qwilcox(0.5-conflevel/2, n1, n2) :
n1*n2/2-qnorm(0.5+conflevel/2)*sqrt(n1*n2*(n1+n2+1)/12)
n = n1*n2
else
if method == "paired.sample"
x = y-x
end
n1 = length(x)
means = (x .+ x') ./ 2
means2 = [means[i, j] for i = 1:n1, j = 1:n1 if j >= i]
vector = sort!(means2)
k = n1 <= 300 ?
qsignrank(0.5-conflevel/2, n1) :
n1*(n1+1)/4-qnorm(0.5+conflevel/2)*sqrt(n1*(n1+1)*(2*n1+1)/24)
n = n1*(n1+1) ÷ 2
end
k = Int(k)
LCL, UCL = vector[[k, n+1-k]]
println("$(Int(100conflevel))% 信頼区間 [$LCL, $UCL]")
Dict(:method => method, :LCL => LCL, :UCL => UCL)
end
# 独立二標本（各 10 ケースずつ）のデータ（ファイルから読んでも良い）
x = [38, 26, 29, 41, 36, 31, 32, 30, 35, 33]
y = [45, 28, 27, 38, 40, 42, 39, 39, 34, 45]
ci4median(x, y; method="independent.sample") # 95% 信頼区間 [-10, 1]
# 対応のある二標本（各 10 ケースずつ）のデータ（ファイルから読んでも良い）
a = [10.6, 5.2, 8.4, 9.0, 6.6, 4.6, 14.1, 5.2, 4.4, 17.4, 7.2]
b = [14.6, 15.6, 20.2, 20.9, 24.0, 25.0, 35.2, 30.2, 30.0, 46.2, 37.0]
ci4median(a, b, method="paired.sample") # 95% 信頼区間 [11.899999999999999, 25.1]
# 一標本のデータ（ファイルから読んでも良い）
ci4median(b-a, method="one.sample") # 95% 信頼区間 [11.899999999999999, 25.1]
ci4median(b-a) # 95% 信頼区間 [11.899999999999999, 25.1]
