#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。
散布図，確率楕円，回帰直線，信頼限界帯，MA regression，RMA regression
http://aoki2.si.gunma-u.ac.jp/R/scatter.html
ファイル名: scatterplot.jl 関数名: scatterplot
翻訳するときに書いたメモ
あれこれオプションを追加すると，引数も増える。
==========#
using Statistics, Distributions, LinearAlgebra, Plots
function scatterplot(x, y;
el =false, el2 =(:black, :solid, 0.5),
lrl=false, lrl2=(:black, :solid, 0.5),
cb =false, cb2 =(:blue, :dot, 0.5),
cb3 =(:brown, :dash, 0.5),
ma =false, ma2 =(:magenta, :solid, 0.5),
rma=false, rma2=(:gray, :solid, 0.5),
marker=(:red, :red, 0.25), factor=0.15,
α=0.05, acc=2000, xlab="", ylab="")
function abline(a, b, col, sty, wid)
plot!([minx, maxx], [a + b * minx, a + b * maxx],
linecolor=col, linestyle=sty, linewidth=wid, label="")
end
function ellipsedraw(draw, el2, α, acc)
λ = eigvals(v)
a = sqrt(vxy^2 / ((λ[1] - vx)^2 + vxy^2))
b = (λ[1] - vx) * a / vxy
θ = atan(a/b)
k = sqrt(-2log(α))
l1 = sqrt(λ[2]) * k
l2 = sqrt(λ[1]) * k
x2 = range(-l1, l1, length=acc)
tmp = [tmp0 < 0 ? 0 : tmp0 for tmp0 in 1 .- x2 .^ 2 ./ l1 ^ 2]
y2 = l2 .* sqrt.(tmp)
x2 = vcat(x2, reverse(x2))
y2 = vcat(y2, -reverse(y2))
s0 = sin(θ)
c0 = cos(θ)
xx = c0 .* x2 .+ s0 .* y2 .+ mean(x)
yy = -s0 .* x2 .+ c0 .* y2 .+ mean(y)
if draw
plot!(xx, yy, linecolor=el2[1], linestyle=el2[2], linewidth=el2[3], label="")
end
return minimum(xx), maximum(xx), minimum(yy), maximum(yy)
end
function conflimit(lrl2, cb, cb2, cb3, α)
b = vxy / vx
a = mean(y) - b * mean(x)
abline(a, b, lrl2[1], lrl2[2], lrl2[3])
if cb
sx2 = vx*(n-1)
x1 = range(minx, maxx, length=2000)
y1 = a .+ b .* x1
ta = cquantile(TDist(n - 2), α/2)
Ve = (vy - vxy^2 / vx) * (n - 1) / (n - 2)
temp = ta * sqrt(Ve) * sqrt.(1/n .+ (x1 .- mean(x)) .^ 2 ./ sx2)
y2 = y1 - temp;
plot!(x1, y2, linecolor=cb2[1], linestyle=cb2[2], linewidth=cb2[3], label="")
y2 = y1 + temp;
plot!(x1, y2, linecolor=cb2[1], linestyle=cb2[2], linewidth=cb2[3], label="")
temp = ta * sqrt(Ve) * sqrt.(1 + 1/n .+ (x1 .- mean(x)) .^ 2 ./ sx2);
y2 = y1 - temp;
plot!(x1, y2, linecolor=cb3[1], linestyle=cb3[2], linewidth=cb3[3], label="")
y2 = y1 + temp;
plot!(x1, y2, linecolor=cb3[1], linestyle=cb3[2], linewidth=cb3[3], label="")
end
print("LS(least squares) ")
println("intercept = $(round(a, digits=5)), slope = $(round(b, digits=5))")
end
function MA(ma2)
b = vxy / (eigvals(v)[2] - vy)
a = mean(y) - b * mean(x)
abline(a, b, ma2[1], ma2[2], ma2[3])
print("MA(major axis) ")
println("intercept = $(round(a, digits=5)), slope = $(round(b, digits=5))")
Dict(:intercept => a, :slope => b)
end
function RMA(rma2)
b = sign(vxy) * sqrt(vy / vx)
a = mean(y) - b * mean(x)
abline(a, b, rma2[1], rma2[2], rma2[3])
print("RMA(reduced major axis) ")
println("intercept = $(round(a, digits=5)), slope = $(round(b, digits=5))")
Dict(:intercept => a, :slope => b)
end
function defminmax(minx0, maxx0, x)
margin = (max(maxx0, maximum(x)) - min(minx0, minimum(x))) * factor
return minx0 - margin, maxx0 + margin
end
pyplot()
v = cov(hcat(x, y))
vx = v[1, 1]
vy = v[2, 2]
vxy = v[1, 2]
minx, maxx, miny, maxy = ellipsedraw(false, el2, α, acc)
minx, maxx = defminmax(minx, maxx, x)
miny, maxy = defminmax(miny, maxy, y)
n = length(x)
p = scatter(x, y, xlims=(minx, maxx), ylims=(miny, maxy),
grid=false, tick_direction=:out, color=marker[1],
markerstrokecolor=marker[2], alpha=marker[3],
xlabel=xlab, ylabel=ylab, label="")
!el || ellipsedraw(true, el2, α, acc)
!lrl || conflimit(lrl2, cb, cb2, cb3, α)
!ma || MA(ma2)
!rma || RMA(rma2)
display(p)
end
x = [132, 146, 140, 196, 132, 154, 154, 168, 140, 140, 156, 114, 134, 116, 150, 178, 150, 120, 150, 146];
y = [90, 90, 84, 96, 90, 90, 74, 92, 60, 82, 80, 62, 80, 80, 76, 98, 86, 70, 80, 80];
scatterplot(x, y, el=true, lrl=true, cb=true)
z = randn((1000, 20));
x2 = mean(z[:, 1:12], dims=2);
y2 = mean(z[:, 8:20], dims=2);
scatterplot(x2, y2, el=true)
