Julia に翻訳--215 非線形回帰

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

非線形曲線回帰
http://aoki2.si.gunma-u.ac.jp/Python/nonlinear_fitting.pdf

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

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

==========#

using LinearAlgebra, NamedArrays, Plots

function nonlinearfitting(x, y; inival=[], model="Asymptotic", method="Marquardt", # or "Simplex"
    maxit=1000, epsilon=1e-7)
    #
    residual(p, x, y) = sum((y .- func(p, x)) .^ 2)
    #
    function marquardt(x, y, inival; maxit, epsilon)
        #
        function solinvfit(p0, x, y)
            n = length(x)
            ip = length(p0)
            ff = zeros(ip)
            fj = zeros(n, ip + 1)
            d = zeros(n + ip)
            delta = zeros(ip)
            p = p0
            resid = residual(p, x, y)
            fval = func(p, x, fj=fj)
            fj[:, ip+1] = y - fval
            ff = fj[:, 1:ip]' * fj[:, ip+1]
            fj = vcat(fj, diagm(-1 => fill(sqrt(cterm), ip))[2:end, :])
            n += ip
            sl = 1e-30
            for j = 1:ip
                d[j:end] = fj[j:end, j]
                h = sum(d[j:end] .^ 2)
                if h < sl
                    g = 0
                else
                    g = sqrt(h)
                    f = fj[j, j]
                    f >= 0 && (g = -g)
                    d[j] = f - g
                    h -= f * g
                    for k = j + 1:ip + 1 # k = j+1
                        s = sum(d[j:end] .* fj[j:end, k]) / h
                        fj[j:end, k] -= d[j:end] .* s
                    end
                end
                fj[j, j] = g
            end
            for i = ip:-1:1 # i = 1
                tot = fj[i, ip + 1] - sum(delta[i + 1:end] .* fj[i, i + 1:ip])
                delta[i] = fj[i, i] == 0 ? 0 : tot / fj[i, i]
            end
            delta, ff, resid
        end
        #
        errset = false
        flag = false
        diverg = 0
        n = length(x)
        ip = ip0 = length(inival)
        p0 = inival
        p = zeros(ip)
        p1 = zeros(ip)
        cterm = 100
        for loop = 1:maxit
            delta, ff, res0 = solinvfit(p0, x, y)
            p1 = p0 + delta
            flag = any(abs.(delta ./ p1) .> epsilon)
            check = -sum((ff .+ cterm .* delta) .* delta)
            p = p1
            res1 = residual(p1, x, y)
            res1 == 0 && break
            rchange = res1 - res0
            if rchange <= 0
                diverg = 0
                errset = false
                check = rchange / check
                if -rchange / res1 < epsilon && flag == false
                    return p
                end
                if check > 0.75
                    cterm *= 0.5
                elseif check < 0.25
                    cterm *= 5
                end
                p0 = p1
            else
                diverg += 1
                errset = true
                if diverg > 10
                    flag && error("not converged")
                    p0 = p1
                else
                    cterm *= 5
                end
            end
        end
        error("not converged")
    end
    #
    function simplex(x, y, inival; maxit, epsilon)
        LO = 0.8
        HI = 1.2
        ip = length(inival)
        ip1 = ip + 1
        ip2 = ip + 2
        ip3 = ip + 3
        pa = float.(reshape(repeat(inival, ip3), ip, :))
        for i = 1:ip
            pa[i, i] = inival[i] * rand(1)[1] * (HI - LO) + LO
        end
        res = vcat([residual(pa[:, i], x, y) for i = 1:ip1], 0, 0)
        for loops = 1:maxit
            res0 = res[1:ip1]
            mx = argmax(res0)
            mi = argmin(res0)
            s = sum(pa[:, 1:ip1], dims = 2)
            if res[mx] < epsilon || res[mi] < epsilon ||
               (res[mx] - res[mi]) / res[mi] < epsilon
               return pa[:, mi]
            end
            i = ip2
            pa[:, ip2] = (2 * s - ip2 * pa[:, mx]) / ip
            res[ip2] = residual(pa[:, ip2], x, y)
            if res[ip2] < res[mi]
                pa[:, ip3] = (3 * s - (2 * ip1 + 1) * pa[:, mx]) / ip
                res[ip3] = residual(pa[:, ip3], x, y)
                res[ip3] <= res[ip2] && (i = ip3)
            elseif res[ip2] > res[mx]
                pa[:, ip3] = s / ip1
                res[ip3] = residual(pa[:, ip3], x, y)
                if res[ip3] >= res[mx]
                    for i = 1:ip1
                        if i != mi
                            pa[:, i] = (pa[:, i] + pa[:, mi]) * 0.5
                            res[i] = residual(pa[:, i], x, y)
                        end
                    end
                    i = 0
                else
                    i = ip3
                end
            end
            if i > 0
                pa[:, mx] = pa[:, i]
                res[mx] = res[i]
            end
        end
        error("not converged!")
    end
    #
    if  model == "Asymptotic"
        func = Asymptotic; ip = 3
        def = "\$y=ab^x+c\$"
    elseif model == "Power2"
        func = Power2; ip = 3
        def = "\$y=ax^b+c\$"
    elseif model == "Exponential1"
        func = Exponential1; ip = 2
        def = "\$y=ab^x\$"
    elseif model == "Power1"
        func = Power1; ip = 2
        def = "\$y=ax^b\$"
    elseif model == "Exponential2"
        func = Exponential2; ip = 2
        def = "\$y=a(1-\\exp(-bx))\$"
    elseif model == "Logistic1"
        func = Logistic1; ip = 3
        def = "\$y=a/(1+b\\exp(-cx))\$"
    elseif model == "Logistic2"
        func = Logistic2; ip = 4
        def = "\$y=a/(1+b\\exp(-cx))+d\$"
    elseif model == "Logistic3"
        func = Logistic3; ip = 3
        def = "\$y=1/(1+a\\exp(-bx))\$"
    elseif model == "Logistic4"
        func = Logistic4; ip = 4
        def = "\$y=a+\\frac{b-a}{1+(c/x)^d}\$"
    elseif model == "Gomperts"
        func = Gomperts; ip = 3
        def = "\$y=ab^\\exp(-cx)\$"
    elseif model == "Weibull"
        func = Weibull; ip = 2
        def = "\$y=1-\\exp(-x^a/b)\$"
    elseif model == "Sin"
        func = Sin; ip = 5
        def = "\$y=a\\sin(bx+c)+dx+e\$"
    elseif model == "Hyperbola1"
        func = Hyperbola1; ip = 3
        def = "\$y=a+\\frac{b}{x+c}\$"
    elseif model == "Hyperbola2"
        func = Hyperbola2; ip = 4
        def = "\$y=\\frac{a}{bx^2+cx+d}\$"
    elseif callable(model)
        func = model
        model = "user defined function"
        def = "\$y\$"
    else
        error("Not yet inplemented")
    end
    length(x) != length(y) && error("lengths of x && y are not equal")
    if model == "user defined function"
        ip = length(inival)
    elseif length(inival) == 0
        inival = fill(1, ip)
    elseif length(inival) != ip
        error("number of initial value must be $ip")
    end
    if method == "Marquardt"
        p = marquardt(x, y, inival; maxit, epsilon)
    else
        p = simplex(x, y, inival; maxit, epsilon)
    end
    predicted = func(p, x)
    resid = sum((predicted .- y) .^ 2)
    println("$model funciton by $(titlecase(method)) method")
    println("estimated parameters $def\n", NamedArray(hcat(Char.(96 .+ (1:ip)), p), (1:ip, ["parameter", "estimated"])))
    println("residual sum of squares $resid")
    println(NamedArray(hcat(x, y, predicted, y .- predicted), (1:length(x), ["x", "y", "pred.", "resid."])))
    Dict(:p => p, :predicted => predicted, :resid => resid, :x => x, :y => y, :model => model, :method => method, :func => func, :def => def)
end

function plotfittedresult(obj)
    #pyplot()
    plt = scatter(obj[:x], obj[:y], grid=false, tick_direction=:out,
        xlabel="\$x\$", ylabel=obj[:def], title=obj[:model], label="")
    minx, maxx = extrema(obj[:x])
    x2 = range(minx, maxx, length=2001)
    y2 = obj[:func](obj[:p], x2)
    plot!(x2, y2, label="")
end

function Asymptotic(p, x; fj=[])
    if length(fj) != 0
        fj[:, 1] = p[2] .^ x
        fj[:, 2] = p[1] .* x .* p[2] .^ (x .- 1)
        fj[:, 3] .= 1
    end
    p[1] .* p[2] .^ x .+ p[3]
end

x = (0:10) / 2 # length(x)
y = [52, 53, 56, 60, 68, 81, 104, 144, 212, 331, 536];
inival=[1,1,1];
obj = nonlinearfitting(x, y, model="Asymptotic")
#=====
Asymptotic funciton by Marquardt method
estimated parameters $y=ab^x+c$
3×2 Named Matrix{Any}
A ╲ B │ parameter  estimated
──────┼─────────────────────
1     │       'a'    2.02823
2     │       'b'    2.99196
3     │       'c'    49.7351
residual sum of squares 0.4215805168503011
11×4 Named Matrix{Float64}
A ╲ B │          x           y       pred.      resid.
──────┼───────────────────────────────────────────────
1     │        0.0        52.0     51.7633    0.236676
2     │        0.5        53.0     53.2434   -0.243384
3     │        1.0        56.0     55.8035    0.196514
4     │        1.5        60.0     60.2318   -0.231767
5     │        2.0        68.0     67.8915    0.108511
6     │        2.5        81.0     81.1407   -0.140729
7     │        3.0       104.0     104.058  -0.0583134
8     │        3.5       144.0       143.7    0.300491
9     │        4.0       212.0     212.268   -0.268009
10    │        4.5       331.0     330.873    0.127114
11    │        5.0       536.0     536.027  -0.0271035
=====#

plotfittedresult(obj) # savefig("fig1.png")
その他の関数への当てはめ例

function Power2(p, x; fj=[])
    temp1 = x .^ p[2]
    if length(fj) != 0
        fj[:, 1] = temp1
        fj[:, 2] = p[1] .* log.(x) .* temp1
        fj[:, 3] .= 1
    end
    p[1] .* temp1 .+ p[3]
end

x = float.(1:10);
y = [7, 21, 59, 133, 255, 436, 691, 1029, 1463, 2005];
obj = nonlinearfitting(x, y, model = "Power2")
plotfittedresult(obj)

function Exponential1(p, x; fj=[])
    if length(fj) != 0
        fj[:, 1] = p[2] .^ x
        fj[:, 2] = p[1] .* x .* p[2] .^ (x .- 1)
    end
    p[1] .* p[2] .^ x
end

x = 1:10;
y = [6, 18, 54, 162, 486, 1458, 4374, 13102, 39366, 118098];
obj = nonlinearfitting(x, y, model="Exponential1")
plotfittedresult(obj)

function Power1(p, x; fj=[])
    temp = x .^ p[2]
    if length(fj) != 0
        fj[:, 1] = temp
        fj[:, 2] = p[1] .* temp .* log.(x)
    end
    p[1] .* temp
end

 x = 1:10;
 y = [2, 16, 54, 128, 250, 432, 680, 1024, 1458, 2000];
 obj = nonlinearfitting(x, y, model="Power1")
 plotfittedresult(obj)

function Exponential2(p, x; fj=[])
    temp = exp.(-x .* p[2])
    if length(fj) != 0
        fj[:, 1] = 1 .- temp
        fj[:, 2] = p[1] .* temp .* x
    end
    p[1] .* (1 .- temp)
end

x = 1:10;
y = [0.518, 0.902, 1.187, 1.398, 1.530, 1.669, 1.755, 1.819, 1.866, 1.900];
obj = nonlinearfitting(x, y, model="Exponential2")
plotfittedresult(obj)

function Logistic1(p, x; fj=[])
    if length(fj) != 0
        temp1 = exp.(p[3] .* x)
        temp2 = exp.(p[3] .* x .* 2)
        temp3 = 1 ./ (temp2 .+ 2p[2] .* temp1 .+ p[2] ^ 2)
        fj[:, 1] = temp1 ./ (temp1 .+ p[2])
        fj[:, 2] = -p[1] .* temp1 .* temp3
        fj[:, 3] = p[1] * p[2] .* x .* temp1 .* temp3
    end
    p[1] ./ (1 .+ p[2] .* exp.(-p[3] .* x))
end

x = 1:10;
y = [1.538, 1.573, 1.606, 1.636, 1.675, 1.692, 1.717, 1.741, 1.762, 1.783];
obj = nonlinearfitting(x, y, model="Logistic1")
plotfittedresult(obj)

function Logistic2(p, x; fj=[])
    if length(fj) != 0
        temp1 = exp.(-p[3] .* x)
        temp2 = 1 .+ p[2] * temp1
        temp3 = 1 ./ temp2 .^ 2
        fj[:, 1] = 1 ./ temp2
        fj[:, 2] = -p[1] .* temp1 .* temp3
        fj[:, 3] = p[1] * p[2] .* temp1 .* x .* temp3
        fj[:, 4] .= 1
    end
    p[1] ./ (1 .+ p[2] * exp.(-p[3] .* x)) .+ p[4]
end

x = 0:10;
y = [2.538, 2.573, 2.606, 2.636, 2.655, 2.692, 2.717, 2.741, 2.762, 2.783, 2.801];
obj = nonlinearfitting(x, y, model="Logistic2")
plotfittedresult(obj)

function Logistic3(p, x; fj=[])
    temp1 = 1 ./ (1 .+ p[1] .* exp.(-p[2] .* x))
    if length(fj) != 0
        temp2 = exp.(-p[2] .* x) .* temp1 .^ 2
        fj[:, 1] = -temp2
        fj[:, 2] = p[1] .* x .* temp2
    end
    temp1
end

x = 0:10;
y = [0.333, 0.403, 0.477, 0.552, 0.604, 0.691, 0.752, 0.803, 0.846, 0.882, 0.909];
obj = nonlinearfitting(x, y, model="Logistic3")
plotfittedresult(obj)

function Logistic4(p, x; fj=[])
    temp1 = p[3] ./ x
    temp = 1 .+ temp1 .^ p[4]
    temp2 = temp .* temp
    if length(fj) != 0
        fj[:, 1] = 1 .- 1 ./ temp
        fj[:, 2] = 1 ./ temp
        fj[:, 3] = (p[1] - p[2]) * p[4] .* temp1 .^ (p[4] - 1) ./ temp2 ./ x
        fj[:, 4] = (p[1] - p[2]) .* temp1 .^ p[4] .* log.(temp1) ./ temp2
    end
    p[1] .+ (p[2] - p[1]) ./ temp
end

x = 1:10;
y = [0.976, 0.816, 0.733, 0.679, 0.631, 0.612, 0.59, 0.571, 0.555, 0.542];
obj = nonlinearfitting(x, y, model="Logistic4")
plotfittedresult(obj)

function Gomperts(p, x; fj=[])
    if length(fj) != 0
        temp1 = exp.(p[3] .* x)
        temp2 = p[2] .^ (1 ./ temp1)
        fj[:, 1] = temp2
        fj[:, 2] = p[1] .* temp2 ./ (p[2] .* temp1)
        fj[:, 3] = -log(p[2]) * p[1] .* x .* temp2 ./ temp1
    end
    p[1] .* p[2] .^ exp.(-p[3] .* x)
end

x = 1:10
y = [0.964, 1.284, 1.529, 1.699, 1.812, 1.884, 1.929, 1.956, 1.973, 1.984];
obj = nonlinearfitting(x, y, model="Gomperts")
plotfittedresult(obj)

function Weibull(p, x; fj=[])
    temp2 = x .^ p[2]
    temp = exp.(-temp2 ./ p[1])
    if length(fj) != 0
        fj[:, 1] = -temp .* temp2 ./ p[1] ^ 2
        fj[:, 2] = temp .* temp2 .* log.(x) ./ p[1]
    end
    1 .- temp
end

x = 1:10;
y = [0.283, 0.418, 0.513, 0.565, 0.642, 0.689, 0.728, 0.76, 0.788, 0.812];
obj = nonlinearfitting(x, y, model="Weibull")
plotfittedresult(obj)

function Sin(p, x; fj=[])
    temp = sin.(p[2] .* x .+ p[3])
    if length(fj) != 0
        temp2 = cos.(p[2] .* x .+ p[3])
        fj[:, 1] = temp
        fj[:, 2] = p[1] .* x .* temp2
        fj[:, 3] = p[1] .* temp2
        fj[:, 4] = x
        fj[:, 5] .= 1
    end
    p[1] .* temp .+ p[4] .* x .+ p[5]
end

x = 1:30;
y = [5.896, 6.539, 5.847, 4.427, 3.23, 3.1, 4.382, 6.754, 9.382,
    11.314, 11.922, 11.203, 9.777, 8.597, 8.5, 9.814, 12.203, 14.826,
    16.731, 17.305, 16.559, 15.127, 13.965, 13.901, 15.247, 17.653,
    20.269, 22.147, 22.686, 21.915];
obj = nonlinearfitting(x, y, model="Sin")
plotfittedresult(obj)

function Hyperbola1(p, x; fj=[])
    temp = x .+ p[3]
    if length(fj) != 0
        fj[:, 1] .= 1
        fj[:, 2] = 1 ./ temp
        fj[:, 3] = -p[2] ./ temp .^ 2
    end
    p[1] .+ p[2] ./ temp
end

x = 1:10;
y = [23.2, 16.5, 11.571, 9.667, 8.455, 7.615, 7, 6.529, 6.158, 5.857];
obj = nonlinearfitting(x, y, model="Hyperbola1")
plotfittedresult(obj)

function Hyperbola2(p, x; fj=[])
    temp = p[2] .* x .^ 2 .+ p[3] .* x .+ p[4]
    if length(fj) != 0
        fj[:, 1] = 1 ./ temp
        temp2 = -p[1] ./ temp .^ 2
        fj[:, 2] = temp2 .* x .^ 2
        fj[:, 3] = temp2 .* x
        fj[:, 4] = temp2
    end
    p[1] ./ temp
end

x = 1:10;
y = [3.922, 1.282, 0.623, 0.396, 0.241, 0.17, 0.127, 0.098, 0.078, 0.063];
obj = nonlinearfitting(x, y, model="Hyperbola2")
plotfittedresult(obj)