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