Julia：使用 CurveFit 包非线性

Question

为什么下面的代码不起作用？

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}) 这可能源于 调用构造函数 Float64(...)，因为类型构造函数下降 返回转换方法。最接近的候选者是：call{T}(::Type{T}, ::Any) convert(::Type{Float64}, !Matched::Int8) 转换(::Type{Float64}, !Matched::Int16)

在加载 In[269] 时，从第 8 行开始的表达式 在 /home/jmarcellopereira/.julia/v0.4/CurveFit/src/nonlinfit.jl:75 的非线性拟合中

Answer 1

fun 函数必须返回 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 函数必须返回单个 Float64，因此运算符不能是“点版本”(.*)。

通过调用nonlinear_fit函数，第一个参数必须是一个数组Nx2，第一列包含x的N个值，第二列包含y的N个值，所以你必须在两列数组中连接两个向量 x 和 y：

xy = [x y]

最后，调用函数：

coefs, converged, iter = CurveFit.nonlinear_fit(xy , fun , a0 , eps, maxiter ) 

回答您关于返回系数的评论不正确：

y = 1 * sin (x * a2-a3) 是一个谐波函数，因此从函数调用返回的系数在很大程度上取决于参数a0（“初步猜测对于每个拟合参数“）您将作为第三个参数发送（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 测试}

Answer 2

来自文档：

我们正在努力适应fun(x, a) = 0的关系

因此，如果您想以以下方式查找a 的元素：对于[x0 y0] 中的每个xi,yi => a[1].*sin(a[2].*xi - a[3])==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);

Answer 3

我发现LsqFit package 使用起来更简单，只需先定义模型并用您的数据“拟合”它：

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 强调的初始条件的依赖性问题..