Python曲线拟合变化点

Question

由于我真的很难从 R 代码转换为 Python 代码，因此我想寻求一些帮助。我想使用的代码是从stackexchange的数学论坛提供给我的。

https://math.stackexchange.com/questions/2205573/curve-fitting-on-dataset

我明白发生了什么。但是我真的很难尝试解决 R 代码，因为我从未见过它。我已经编写了返回平方和的函数。但我被困在如何使用类似于 optim 函数的函数。而且我真的不喜欢对初始值的猜测。我希望更好地运行并重新运行一种优化函数，直到我得到想要的结果，因为我对近乎完美的曲线拟合的需求非常高。

def model (par,x):
    n = len(x)
    res = []
    for i in range(1,n):
        A0 = par[3] + (par[4]-par[1])*par[6] + (par[5]-par[2])*par[6]**2
        if(x[i] == par[6]):
            res[i] = A0 + par[1]*x[i] + par[2]*x[i]**2
        else:
            res[i] = par[3] + par[4]*x[i] + par[5]*x[i]**2
    return res

这是我的模型函数...

def sum_squares (par, x, y):
    ss = sum((y-model(par,x))^2)
    return ss

这是平方和

但我不知道如何转换：

 #I found these initial values with a few minutes of guess and check.
 par0 <- c(7,-1,-395,70,-2.3,10)
 sol <- optim(par= par0, fn=sqerror, x=x, y=y)$par

到 Python 代码...

Answer 1

我编写了一个开源 Python 包（BSD 许可证），它具有 scipy Levenberg-Marquardt 求解器的遗传算法（差分进化）前端，它的功能与您在问题中描述的类似。 github地址是：

https://github.com/zunzun/pyeq3

它带有一个相当容易使用的“用户定义函数”示例：

https://github.com/zunzun/pyeq3/blob/master/Examples/Simple/FitUserDefinedFunction_2D.py

以及命令行、GUI、集群、并行和基于 Web 的示例。您可以使用“pip3 install pyeq3”安装该软件包，看看它是否适合您的需求。

Answer 2

看来我已经能够解决问题了。

def model (par,x):
    n = len(x)
    res = np.array([]) 
    for i in range(0,n):
        A0 = par[2] + (par[3]-par[0])*par[5] + (par[4]-par[1])*par[5]**2
        if(x[i] <= par[5]):
            res = np.append(res, A0 + par[0]*x[i] + par[1]*x[i]**2)
        else:
            res = np.append(res,par[2] + par[3]*x[i] + par[4]*x[i]**2)
    return res

def sum_squares (par, x, y):
    ss = sum((y-model(par,x))**2)
    print('Sum of squares = {0}'.format(ss))
    return ss

然后我使用了如下函数：

parameter = sy.array([0.0,-8.0,0.0018,0.0018,0,200])
res = least_squares(sum_squares, parameter, bounds=(-360,360), args=(x1,y1),verbose = 1)

唯一的问题是它不会产生我正在寻找的结果......这主要是因为我的 x 值是 [0,360] 并且 Y 值仅变化大约 0.2，所以这是一个难题破解这个函数，它会产生这个（糟糕的）结果：

Result

Answer 3

我认为 x 值 [0, 360] 和 y 值（你说是 ~0.2）的范围可能不是问题。获得良好的参数初始值可能更为重要。

在带有 numpy / scipy 的 Python 中，您肯定希望 不 循环 x 的值，而是做一些更像

def model(par,x):
    res = par[2] + par[3]*x + par[4]*x**2        
    A0  = par[2] + (par[3]-par[0])*par[5] + (par[4]-par[1])*par[5]**2
    res[np.where(x <= par[5])] = A0 + par[0]*x + par[1]*x**2 
    return res

我不清楚这种形式是否真的是你想要的：为什么 A0（一个独立于 x 的值添加到模型的一部分）如此复杂并且相互依赖于其他参数？

更重要的是，你的sum_of_squares() 函数实际上不是least_squares() 想要的：你应该返回残差数组，你不应该自己做平方和。所以，应该是

def sum_of_squares(par, x, y): 
    return (y - model(par, x))

但最重要的是，有一个概念问题可能会困扰此模型：您的 par[5] 旨在表示模型更改形式的断点。对于这些优化例程来说，这将是非常难以找到的。这些例程通常对每个参数值进行非常小的更改，以估计残差数组相对于该变量的导数，以便弄清楚如何更改该变量。对于本质上用作整数的参数，初始值的微小变化将完全没有影响，算法将无法确定该参数的值。使用一些 scipy.optimize 算法（特别是 leastsq），您可以指定要进行的相对更改的比例。使用leastsq 称为epsfcn。您可能需要将其设置为 0.3 或 1.0 才能使断点正常工作。不幸的是，这不能按变量设置，只能按拟合设置。您可能需要尝试 least_squares 或 leastsq 的此选项和其他选项。