将数据拟合到 post 处理后的非线性模型时,curve_fit 的算法是什么?

What is curve_fit's algorithm when fitting data to a post processed non-linear model?

我是编程工具新手,Python,对数值计算了解甚少。我试图理解 Curve_fit 的源代码,但对我来说太难了。

如果我们直接对非线性模型进行曲线拟合,通过计算Hessian矩阵,我想我们可以找到所有参数的临界点,不管局部还是全局,最小值还是最大值。

但是,如果包含模型的函数无法微分,curve_fit 如何实际计算答案?

我想算法是在所有数据点的每个参数的每个区间中找到 Chi^2 的局部最小值。如果这样做,默认间隔有多大,每次参数试验之间的距离是多少,最大迭代次数是多少?如果像我说的那样,这个迭代在多个参数之间是如何进行的呢?他们每次都作为一组单独尝试吗?

我想了解算法,以便我可以编写一个好的函数来应用curve_fit。我尝试适应的函数中存在许多脏修复,不同的小更改出现不同的错误,因此我无法将代码放在这里,因为我不知道问题所在。

此外,关于 sigma 输入,如果我不知道 y 误差,如果默认 sigma 与数据值相比太大,结果会怎样?或者,如果模型的灵活性对 post 过程影响不大,或者 sigma 远小于数据到函数的距离,会发生什么情况?

要开始了解 scipy curve_fit 如何尝试解决问题,请先阅读非线性优化方法,例如 Levenberg-Marquardt (LM) 方法(请参阅 https://en.wikipedia.org/wiki/Levenberg%E2%80%93Marquardt_algorithm).这是 curve_fit 在未指定参数值范围时使用的方法,它会给出您寻求的大部分答案。

此方法(以及 scipy.optimize 中的其他求解器)从初始值开始并尝试优化这些值以找到最小的结果残差数组(在最小二乘意义上:要清楚,该方法需要一个残差数组并自行进行平方和求和)。 为了找到移动参数值的方向,它使用一阶导数(残差与变量的关系)或雅可比行列式。对此没有解析形式是很常见的,在这种情况下,数值导数是通过在参数值中采取非常小的步长(通常围绕机器精度)来计算的。

有了这些导数和最后的步长,LM 方法决定了要找到最小残差的步长。将残差视为每个参数的二次函数是一种常见的(通常是很好的)近似。如果该值远离底部,则采取导数中的线性步骤(即,仅遵循局部斜率)是非常好的。但在接近最终解决方案时,使用二阶导数(又名 Hessian、曲率或某种意义上的 "history")会很有帮助。 LM 将这两种方法结合起来,试图快速找到解决方案。在每一步,它首先计算导数(和第一步之后的二阶导数),然后计算所有参数的步长。它重复此过程,直到解决方案(卡方)与公差相比没有改善。在 scipy 的 leastsq(通常被 curve_fit 使用)中,可以指定用于计算数值导数的步长和停止拟合的公差。

通常,当使用 curve_fitleastsqleast_squares 时,您无需担心有关如何找到解决方案的任何这些细节。好吧,除了如果您 有解析导数,它可以提高过程的速度和准确性。您的 "model function" 应该只接受传入的参数并计算数据的预测模型。更改传入的值会使算法混淆 - 运行 您的函数具有它知道的值,以便尝试找到解决方案。

其他问题:
如果您不知道数据中的标准错误 (y),请不要太担心。报告卡方的绝对值可能不接近 ndata-nvarys(符合良好拟合的预期),但相对于其他拟合的值应该没问题。

LM 方法的一个特点(至少在 MINPACK -> leastsq -> curve_fit 实现中)是它报告可用于确定不确定性的最终协方差矩阵和拟合参数的相关性。通常报告的值应该将卡方增加 1(即 1-sigma 标准误差)。由于数据的不确定性非常普遍,scipy curve_fit 缩放此协方差以给出值,就好像卡方等于 ndata-nvarys 一样。也就是说,如果您认为拟合良好(并且 sqrt(chi-square/(ndata-nvays)) 可用于估计数据中的误差),参数中的标准误差将非常好。