3

why the code below does not work? 

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

4

3 回答 3

3

该fun函数必须返回r类型为 的残差值Float64，在数据的每次迭代中计算，如下所示：

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，因此运算符不能是“点版本”（.*）。

通过调用非线性函数，第一个参数必须是一个数组 Nx2，第一列包含 N 个值，x第二列包含 N 个值y，因此您必须将两个向量连接起来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 进行测试

于 2016-01-24T07:51:39.917 回答
2

来自文档：

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

因此，如果您想以a以下方式查找元素：对于=>xi,yi中的每个元素，那么正确的方法是：[x0 y0]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);
于 2016-01-18T03:32:18.747 回答
2

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

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

于 2017-08-29T09:59:00.677 回答