#==========
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)
