Python 变化点曲线拟合

Question

因为我真的很难从 R 代码到 Python 代码，所以我想寻求一些帮助。我想使用的代码已经从 stackexchange 的数学论坛提供给我。

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

我明白是怎么回事了。但我真的很难尝试解决 R 代码，因为我从未见过任何东西。我已经将函数编写为 return 平方和。但我坚持如何使用类似于 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

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

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 2

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

https://github.com/zunzun/pyeq3

它附带了一个相当容易使用的 "user-defined function" 示例：

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

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

Answer 3

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

在 Python 和 numpy / scipy 中，您肯定希望 而不是 遍历 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()想要的：你应该return残差数组，你不应该自己做平方和。所以，那应该是

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 的这个选项和其他选项。

Python 变化点曲线拟合

Python curve fit with change point

python

curve

r

curve-fitting