Julia：使用 CurveFit 包非线性
Julia: using CurveFit package non linear
为什么下面的代码不起作用？
xa = [0 0.200000000000000 0.400000000000000 1.00000000000000 1.60000000000000 1.80000000000000 2.00000000000000 2.60000000000000 2.80000000000000 3.00000000000000 3.80000000000000 4.80000000000000 5.00000000000000 5.20000000000000 6.00000000000000 6.20000000000000 7.40000000000000 7.60000000000000 7.80000000000000 8.60000000000000 8.80000000000000 9.00000000000000 9.20000000000000 9.40000000000000 10.0000000000000 10.6000000000000 10.8000000000000 11.2000000000000 11.6000000000000 11.8000000000000 12.2000000000000 12.4000000000000];
ya = [-0.183440428023042 -0.131101157495126 0.0268875670852843 0.300000000120000 0.579048247883555 0.852605831272159 0.935180993484717 1.13328608090532 1.26893326843583 1.10202945535186 1.09201137189664 1.14279083803453 0.811302535321072 0.909735376251797 0.417067545528244 0.460107770989798 -0.516307074859654 -0.333994077331822 -0.504124744955962 -0.945794320817293 -0.915934553082780 -0.975458595671737 -1.09943707404275 -1.11254211607374 -1.29739980589100 -1.23440439602665 -0.953807504156356 -1.12240274852172 -0.609284630192522 -0.592560286759450 -0.402521296049042 -0.510090363150962];
x0 = vec(xa)
y0 = vec(ya)
fun(x,a) = a[1].*sin(a[2].*x - a[3])
a0 = [1,2,3]
eps = 0.000001
maxiter=200
coefs, converged, iter = CurveFit.nonlinear_fit(x0 , fun , a0 , eps, maxiter )
y0b = fit(x0)
Winston.plot(x0, y0, "ob", x0, y0b, "r-", linewidth=3)
Error: LoadError: MethodError: convert has no method matching convert(::Type{Float64}, ::Array{Float64,1}) This may have arisen from
a call to the constructor Float64(...), since type constructors fall
back to convert methods. Closest candidates are: call{T}(::Type{T},
::Any) convert(::Type{Float64}, !Matched::Int8)
convert(::Type{Float64}, !Matched::Int16)
while loading In[269], in expression starting on line 8
in nonlinear_fit at /home/jmarcellopereira/.julia/v0.4/CurveFit/src/nonlinfit.jl:75
来自医生：
we are trying to fit the relationship fun(x, a) = 0
所以，如果你想以这样的方式查找 a 的元素：对于 [x0 y0] => a[1].*sin(a[2].*xi - a[3])==yi 中的每个 xi,yi，那么正确的方法是:
fun(xy,a) = a[1].*sin(a[2].*xy[1] - a[3])-xy[2];
xy=hcat(x0,y0);
coefs,converged,iter = CurveFit.nonlinear_fit(xy,fun,a0,eps,maxiter);
fun函数必须return一个Float64类型的残差值r，在数据的每次迭代中计算，如下：
r = y - fun(x, coefs)
因此您的函数 y=a1*sin(x*a2-a3) 将定义为：
fun(x,a) = x[2]-a[1]*sin(a[2]*x[1] - a[3])
其中：
x[2] is a value of 'y' vector
x[1] is a value of 'x' vector
a[...] is the set of parameters
fun 函数必须 return 单个 Float64，因此运算符不能是 'dot version' (.*).
通过调用nonlinear_fit函数，第一个参数必须是一个数组Nx2，第一列包含x的N个值，第二列，包含 y 的 N 个值，因此您必须在两列数组中连接两个向量 x 和 y：
xy = [x y]
最后，调用函数：
coefs, converged, iter = CurveFit.nonlinear_fit(xy , fun , a0 , eps, maxiter )
回答您关于 returned 系数不正确的评论：
y = 1 * sin (x * a2-a3) 是一个 harmonic 函数，因此函数调用的系数 return 很大程度上取决于 在参数 a0 ("initial guess for each fitting parameter") 上，您将作为第三个参数发送（带有 maxiter=200_000）：
a0=[1.5, 1.5, 1.0]
coefficients: [0.2616335317043578, 1.1471991302529982,0.7048665905560775]
a0=[100.,100.,100.]
coefficients: [-0.4077952060368059, 90.52328921205392, 96.75331155303707]
a0=[1.2, 0.5, 0.5]
coefficients: [1.192007321713507, 0.49426296880933257, 0.19863645732313934]
我认为您得到的结果是谐波，如图所示：
其中：
blue line:
f1(xx)=0.2616335317043578*sin(xx*1.1471991302529982-0.7048665905560775)
yellow line:
f2(xx)=1.192007321713507*sin(xx*0.49426296880933257-0.19863645732313934)
pink line:
f3(xx)=-0.4077952060368059*sin(xx*90.52328921205392-96.75331155303707)
blue dots are your initial data.
图表是用 Gadfly 生成的：
plot(layer(x=x,y=y,Geom.point),layer([f1,f2,f3],0.0, 15.0,Geom.line))
使用 Julia 版本 0.4.3 测试
我发现 LsqFit package 使用起来更简单一些，只需先定义模型，然后 "fit it" 使用您的数据：
using DataFrames, Plots, LsqFit
xa = [0 0.200000000000000 0.400000000000000 1.00000000000000 1.60000000000000 1.80000000000000 2.00000000000000 2.60000000000000 2.80000000000000 3.00000000000000 3.80000000000000 4.80000000000000 5.00000000000000 5.20000000000000 6.00000000000000 6.20000000000000 7.40000000000000 7.60000000000000 7.80000000000000 8.60000000000000 8.80000000000000 9.00000000000000 9.20000000000000 9.40000000000000 10.0000000000000 10.6000000000000 10.8000000000000 11.2000000000000 11.6000000000000 11.8000000000000 12.2000000000000 12.4000000000000];
ya = [-0.183440428023042 -0.131101157495126 0.0268875670852843 0.300000000120000 0.579048247883555 0.852605831272159 0.935180993484717 1.13328608090532 1.26893326843583 1.10202945535186 1.09201137189664 1.14279083803453 0.811302535321072 0.909735376251797 0.417067545528244 0.460107770989798 -0.516307074859654 -0.333994077331822 -0.504124744955962 -0.945794320817293 -0.915934553082780 -0.975458595671737 -1.09943707404275 -1.11254211607374 -1.29739980589100 -1.23440439602665 -0.953807504156356 -1.12240274852172 -0.609284630192522 -0.592560286759450 -0.402521296049042 -0.510090363150962];
x0 = vec(xa)
y0 = vec(ya)
xbase = collect(linspace(minimum(x0),maximum(x0),100))
p0 = [1.2, 0.5, 0.5] # initial value of parameters
fun(x0, p) = p[1] .* sin.(p[2] .* x0 .- p[3]) # definition of the model
fit = curve_fit(fun,x0,y0,p0) # actual fitting job
yFit = [fit.param[1] * sin(fit.param[2] * x - fit.param[3]) for x in xbase] # building the fitted values
# Plotting..
scatter(x0, y0, label="obs")
plot!(xbase, yFit, label="fitted")
请注意，使用 LsqFit 并不能解决 Gomiero 强调的初始条件的依赖问题。
为什么下面的代码不起作用？
xa = [0 0.200000000000000 0.400000000000000 1.00000000000000 1.60000000000000 1.80000000000000 2.00000000000000 2.60000000000000 2.80000000000000 3.00000000000000 3.80000000000000 4.80000000000000 5.00000000000000 5.20000000000000 6.00000000000000 6.20000000000000 7.40000000000000 7.60000000000000 7.80000000000000 8.60000000000000 8.80000000000000 9.00000000000000 9.20000000000000 9.40000000000000 10.0000000000000 10.6000000000000 10.8000000000000 11.2000000000000 11.6000000000000 11.8000000000000 12.2000000000000 12.4000000000000];
ya = [-0.183440428023042 -0.131101157495126 0.0268875670852843 0.300000000120000 0.579048247883555 0.852605831272159 0.935180993484717 1.13328608090532 1.26893326843583 1.10202945535186 1.09201137189664 1.14279083803453 0.811302535321072 0.909735376251797 0.417067545528244 0.460107770989798 -0.516307074859654 -0.333994077331822 -0.504124744955962 -0.945794320817293 -0.915934553082780 -0.975458595671737 -1.09943707404275 -1.11254211607374 -1.29739980589100 -1.23440439602665 -0.953807504156356 -1.12240274852172 -0.609284630192522 -0.592560286759450 -0.402521296049042 -0.510090363150962];
x0 = vec(xa)
y0 = vec(ya)
fun(x,a) = a[1].*sin(a[2].*x - a[3])
a0 = [1,2,3]
eps = 0.000001
maxiter=200
coefs, converged, iter = CurveFit.nonlinear_fit(x0 , fun , a0 , eps, maxiter )
y0b = fit(x0)
Winston.plot(x0, y0, "ob", x0, y0b, "r-", linewidth=3)
Error: LoadError: MethodError:
converthas no method matching convert(::Type{Float64}, ::Array{Float64,1}) This may have arisen from a call to the constructor Float64(...), since type constructors fall back to convert methods. Closest candidates are: call{T}(::Type{T}, ::Any) convert(::Type{Float64}, !Matched::Int8) convert(::Type{Float64}, !Matched::Int16)while loading In[269], in expression starting on line 8 in nonlinear_fit at /home/jmarcellopereira/.julia/v0.4/CurveFit/src/nonlinfit.jl:75
来自医生：
we are trying to fit the relationship
fun(x, a) = 0
所以，如果你想以这样的方式查找 a 的元素：对于 [x0 y0] => a[1].*sin(a[2].*xi - a[3])==yi 中的每个 xi,yi，那么正确的方法是:
fun(xy,a) = a[1].*sin(a[2].*xy[1] - a[3])-xy[2];
xy=hcat(x0,y0);
coefs,converged,iter = CurveFit.nonlinear_fit(xy,fun,a0,eps,maxiter);
fun函数必须return一个Float64类型的残差值r，在数据的每次迭代中计算，如下：
r = y - fun(x, coefs)
因此您的函数 y=a1*sin(x*a2-a3) 将定义为：
fun(x,a) = x[2]-a[1]*sin(a[2]*x[1] - a[3])
其中：
x[2] is a value of 'y' vector
x[1] is a value of 'x' vector
a[...] is the set of parameters
fun 函数必须 return 单个 Float64，因此运算符不能是 'dot version' (.*).
通过调用nonlinear_fit函数，第一个参数必须是一个数组Nx2，第一列包含x的N个值，第二列，包含 y 的 N 个值，因此您必须在两列数组中连接两个向量 x 和 y：
xy = [x y]
最后，调用函数：
coefs, converged, iter = CurveFit.nonlinear_fit(xy , fun , a0 , eps, maxiter )
回答您关于 returned 系数不正确的评论：
y = 1 * sin (x * a2-a3) 是一个 harmonic 函数，因此函数调用的系数 return 很大程度上取决于 在参数 a0 ("initial guess for each fitting parameter") 上，您将作为第三个参数发送（带有 maxiter=200_000）：
a0=[1.5, 1.5, 1.0]
coefficients: [0.2616335317043578, 1.1471991302529982,0.7048665905560775]
a0=[100.,100.,100.]
coefficients: [-0.4077952060368059, 90.52328921205392, 96.75331155303707]
a0=[1.2, 0.5, 0.5]
coefficients: [1.192007321713507, 0.49426296880933257, 0.19863645732313934]
我认为您得到的结果是谐波，如图所示：
其中：
blue line:
f1(xx)=0.2616335317043578*sin(xx*1.1471991302529982-0.7048665905560775)
yellow line:
f2(xx)=1.192007321713507*sin(xx*0.49426296880933257-0.19863645732313934)
pink line:
f3(xx)=-0.4077952060368059*sin(xx*90.52328921205392-96.75331155303707)
blue dots are your initial data.
图表是用 Gadfly 生成的：
plot(layer(x=x,y=y,Geom.point),layer([f1,f2,f3],0.0, 15.0,Geom.line))
使用 Julia 版本 0.4.3 测试
我发现 LsqFit package 使用起来更简单一些，只需先定义模型，然后 "fit it" 使用您的数据：
using DataFrames, Plots, LsqFit
xa = [0 0.200000000000000 0.400000000000000 1.00000000000000 1.60000000000000 1.80000000000000 2.00000000000000 2.60000000000000 2.80000000000000 3.00000000000000 3.80000000000000 4.80000000000000 5.00000000000000 5.20000000000000 6.00000000000000 6.20000000000000 7.40000000000000 7.60000000000000 7.80000000000000 8.60000000000000 8.80000000000000 9.00000000000000 9.20000000000000 9.40000000000000 10.0000000000000 10.6000000000000 10.8000000000000 11.2000000000000 11.6000000000000 11.8000000000000 12.2000000000000 12.4000000000000];
ya = [-0.183440428023042 -0.131101157495126 0.0268875670852843 0.300000000120000 0.579048247883555 0.852605831272159 0.935180993484717 1.13328608090532 1.26893326843583 1.10202945535186 1.09201137189664 1.14279083803453 0.811302535321072 0.909735376251797 0.417067545528244 0.460107770989798 -0.516307074859654 -0.333994077331822 -0.504124744955962 -0.945794320817293 -0.915934553082780 -0.975458595671737 -1.09943707404275 -1.11254211607374 -1.29739980589100 -1.23440439602665 -0.953807504156356 -1.12240274852172 -0.609284630192522 -0.592560286759450 -0.402521296049042 -0.510090363150962];
x0 = vec(xa)
y0 = vec(ya)
xbase = collect(linspace(minimum(x0),maximum(x0),100))
p0 = [1.2, 0.5, 0.5] # initial value of parameters
fun(x0, p) = p[1] .* sin.(p[2] .* x0 .- p[3]) # definition of the model
fit = curve_fit(fun,x0,y0,p0) # actual fitting job
yFit = [fit.param[1] * sin(fit.param[2] * x - fit.param[3]) for x in xbase] # building the fitted values
# Plotting..
scatter(x0, y0, label="obs")
plot!(xbase, yFit, label="fitted")
请注意，使用 LsqFit 并不能解决 Gomiero 强调的初始条件的依赖问题。