Julia に翻訳--044 Reduced Major Axis regression，RMA

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

Reduced Major Axis regression
http://aoki2.si.gunma-u.ac.jp/R/RMA.html

ファイル名: rma.jl　　関数名: rma

翻訳するときに書いたメモ

もはや，特になし。

==========#

using Statistics, Distributions, Printf, Plots

function rma(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
    n = length(x)
    n1 = n-1
    df = n-2
    slope = sign(cor(x, y))*std(y)/std(x)
    intercept = mean(y)-slope*mean(x)
    MSE = (var(y)-cov(x, y)^2/var(x))*n1/df
    SEintercept = sqrt(MSE*(1/n+mean(x)^2/var(x)/n1))
    SEslope = sqrt(MSE/var(x)/n1)
    alpha = (1-sig)/2
    df = n - 2
    t = quantile(TDist(df), [alpha, 1-alpha])
    CLintercept = intercept .+ t .* SEintercept
    CLslope = slope .+ t .* SEslope
    @printf("            Estimate      S.E.     %6.3f%%   %6.3f%%\n", alpha, 1-alpha)
    @printf("Intercept:  %8.7g  %8.7g  [%8.7g, %8.7g]\n", intercept, SEintercept, CLintercept[1], CLintercept[2])
    @printf("Slope:      %8.7g  %8.7g  [%8.7g, %8.7g]\n", slope, SEslope, CLslope[1], CLslope[2])
    plotrma(x, y, intercept, slope, col1, col2, xlab, ylab, legend)
    Dict(:intercept => intercept, :SEintercept => SEintercept,
         :CLintercept => CLintercept, :slope => slope,
         :SEslope => SEslope, :CLslope => CLslope)
end

function plotrma(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, "Reduced Major Axis")
    a, b = lm(x, y)
    abline(a, b, col2, "linear regression")
    display(plt)
end

y = [61, 37, 65, 69, 54, 93, 87, 89, 100, 90, 97]
x = [14, 17, 24, 25, 27, 33, 34, 37, 40, 41, 42]

a = rma(x, y, legend=:top)

            Estimate      S.E.      0.025%    0.975%
Intercept:  12.19378  10.54975  [-11.67141, 36.05898]
Slope:      2.119366  0.3324959  [1.367208, 2.871524]

a = rma(x, y, sig=0.9, legend=:top)

            Estimate      S.E.      0.050%    0.950%
Intercept:  12.19378  10.54975  [-7.145098, 31.53267]
Slope:      2.119366  0.3324959  [1.509864, 2.728869]

Julia に翻訳--043 二本の直線による折れ線回帰

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

二本の直線による折れ線回帰
http://aoki2.si.gunma-u.ac.jp/R/oresen.html

ファイル名: oresen.jl　　関数名: oresen

翻訳するときに書いたメモ

データにより FUNC に何を選ぶかにより違いが出てくるのだろう。

==========#

using Statistics, Printf, Plots, Optim

function oresen(x, y; xlab="", ylab="", col1=:red, col2=:blue,
    FUNC=NelderMead) #FUNC=NelderMead, LBFGS or BFGS =#)
    function ss(par)
        a, b, c, d = par
        flag = x .< a
        xl = x[flag]
        xr = x[.!flag]
        (length(xl) != 0 && length(xr) != 0) || return Inf
        yl = y[flag]
        yr = y[.!flag]
        yle = c .* (xl .- a) .+ b
        yre = d .* (xr .- a) .+ b
        sum((yl .- yle) .^ 2) + sum((yr .- yre) .^ 2)
    end

    function printoresen(par)
        @printf("残差平方和 = %.7g\n", ss(par))
        @printf("交点座標 = ( %.7g, %.7g )\n", par[1], par[2])
        @printf("切片 = %.7g, 傾き = %.7g\n", -par[3]*par[1]+par[2], par[3])
        @printf("切片 = %.7g, 傾き = %.7g\n", -par[4]*par[1]+par[2], par[4])
    end

    function plotoresen(x, y, par, col1, col2)
        pyplot()
        yhat(slope, x) = slope*(x - par[1]) + par[2]
        flag = x .< par[1]
        scatter(x[flag], y[flag], grid=false, tick_direction=:outer,
             xlabel=xlab, ylabel=ylab, color=col1, label="")
        scatter!(x[.!flag], y[.!flag], color=col2, label="")
        minx, maxx = extrema(x)
        margin = (maxx - minx) * 0.05
        minx -= margin
        maxx += margin
        plot!([minx, maxx], [yhat(par[3], minx), yhat(par[3], maxx)],
                color=col1, label="")
        plot!([minx, maxx], [yhat(par[4], minx), yhat(par[4], maxx)],
                color=col2, label="")
    end
    par = [mean(x), mean(y), 1, 1]
    result = optimize(ss, par, FUNC())
    par = Optim.minimizer(result)
    printoresen(par)
    plotoresen(x, y, par, col1, col2)
end

x = [1, 2, 3, 6, 6.1, 4, 5, 7, 8, 9, 10];
y = [2.3, 4.1, 6.4, 12, 12.3, 7.8, 9.7, 14.4, 17.7, 21.5, 24.1];
oresen(x, y)
残差平方和 = 0.7303773
交点座標 = ( 5.85193, 11.33613 )
切片 = 0.5099524, 傾き = 1.850019
切片 = -6.761783, 傾き = 3.09264

oresen(x, y, FUNC=LBFGS)
残差平方和 = 0.5425213
交点座標 = ( 6.528412, 12.93847 )
切片 = 0.312848, 傾き = 1.933951
切片 = -8.54, 傾き = 3.29

Julia の小ネタ--009 最適化，Optimization

Optimization は

# https://julianlsolvers.github.io/Optim.jl/stable/
using Optim
rosenbrock(x) =  (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
result = optimize(rosenbrock, zeros(2), BFGS())
summary(result)
Optim.minimizer(result)
Optim.minimum(result)
Optim.converged(result)

最初 optimize の結果に推定値が出ていないので焦った。

マニュアルを読み進めて，Potim.minimizer() で出力できることがわかった。

Julia に翻訳--042 定点を通る直線回帰式の傾き

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

定点を通る直線回帰式の傾き
http://aoki2.si.gunma-u.ac.jp/R/sregc.html

ファイル名: sregc.jl　　関数名: sregc

翻訳するときに書いたメモ

cx，cy の記述法を変える。

==========#

using Statistics

function sregc(x, y; cx=mean(x), cy=mean(y))
    n = length(x)
    sxy = sum(x .* y)
    sxx = sum(x .^ 2)
    sx = sum(x)
    sy = sum(y)
    slope = (sxy-cy*sx-cx*sy+n*cx*cy)/(sxx-2*cx*sx+n*cx^2)
    intercept = cy - slope * cx
    return (intercept, slope)
end

x = [1,2,3,4,5,6]
y = [3,2,4,1,4,3]

sregc(x, y, cx=0, cy=2) # (2.0, 0.2087912087912088)
sregc(x, y, cx=0, cy=0) # (0.0, 0.6703296703296703)
sregc(x, y) # (2.5333333333333337, 0.08571428571428572)