#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。
Major Axis regression（主成分回帰）
http://aoki2.si.gunma-u.ac.jp/R/MA.html
ファイル名: ma.jl 関数名: ma
翻訳するときに書いたメモ
LinearAlgebra の固有値は降順ではなく昇順
==========#
using Statistics, Distributions, LinearAlgebra, Printf, Plots
function ma(x, y; sig=0.95, col1=:black, col2=:red, xlab="", ylab="",
legend=:best)
# legend :right, :left, :top, :bottom, :inside, :best, :legend,
# :topright, :topleft, :bottomleft, :bottomright
s2 = cov(hcat(x, y))
ev = eigvals(s2)
slope = s2[1, 2]/(ev[2]-s2[2, 2])
intercept = mean(y)-slope*mean(x)
n = length(x)
H = quantile(FDist(1, n-2), sig)/(ev[2]/ev[1]+ev[1]/ev[2]-2)/(n-2)
CL = tan.(atan(slope) .+ [-0.5, 0.5] .* asin(2*sqrt(H)))
alpha = (1-sig)/2
@printf(" Estimate %6.3f%% %6.3f%%\n", alpha, 1-alpha)
@printf("Intercept: %8.7g --- ---\n", intercept)
@printf("Slope: %8.7g [%8.7g, %8.7g]\n", slope, CL[1], CL[2])
plotma(x, y, intercept, slope, col1, col2, xlab, ylab, legend)
Dict(:intercept => intercept, :slope=>slope, :lcl => CL[1], :ucl => CL[2])
end
function plotma(x, y, intercept, slope, col1, col2, xlab, ylab, legend)
abline(intercept, slope, col, label) =
plot!([minx, maxx],
[slope * minx + intercept, slope * maxx + intercept],
color=col, label=label, legend=legend)
function lm(x, y)
b = cov(x, y)/var(x)
return mean(y) - b * mean(x), b
end
minx, maxx = extrema(x)
margin = (maxx - minx) * 0.05
minx -= margin
maxx += margin
pyplot()
plt = scatter(x, y, xlabel=xlab, ylabel=ylab, color=col1, label="")
abline(intercept, slope, col1, "Major Axis regression")
a, b = lm(x, y)
abline(a, b, col2, "linear regression")
display(plt)
end
x = [14.40, 15.20, 11.30, 2.50, 22.70, 14.90, 1.41, 15.81, 4.19, 15.39, 17.25, 9.52]
y = [159, 179, 100, 45, 384, 230, 100, 320, 80, 220, 320, 210]
ma(x, y)
Estimate 0.025% 0.975%
Intercept: -32.55209 --- ---
Slope: 18.93633 [13.43925, 31.99261]
