在 Python 中拟合参数曲线

Posted 2023-03-12

【中文标题】在 Python 中拟合参数曲线

【英文标题】：Fitting Parametric Curves in Python

【发布时间】：2015-11-15 00:14:06

【问题描述】：

我有(X,Y) 形式的实验数据和(x(t;*params),y(t;*params)) 形式的理论模型，其中t 是一个物理（但不可观察）变量，*params 是我想要确定的参数。 t是一个连续变量，模型中x和t之间、y和t之间存在1:1的关系。

在一个完美的世界中，我会知道T 的值（参数的实际值），并且能够进行非常基本的最小二乘拟合来找到*params 的值。 （请注意，我不是试图在我的情节中“连接”x 和 y 的值，就像在 31243002 或 31464345 中一样。）我不能保证在我的真实数据，潜在值T 是单调的，因为我的数据是跨多个周期收集的。

我对手动进行曲线拟合不是很有经验，并且必须使用极其粗糙的方法，而无法轻松访问基本的 scipy 函数。我的基本方法包括：

选择*params 的某个值并将其应用于模型 取t 值的数组并将其放入模型中以创建model(*params) = (x(*params),y(*params)) 数组 将X（数据值）插入model，得到Y_predicted 在Y 和Y_predicted 之间进行最小二乘（或其他）比较 再次为一组新的*params 最终，为*params 选择最佳值

这种方法有几个明显的问题。

1) 我没有足够的编码经验来开发一个非常好的“再做一次”，而不是“尝试解决方案空间中的所有内容”，也许是“在粗略的网格中尝试所有内容”，然后“再次尝试所有内容”在粗网格的热点中稍细的网格中。”我尝试使用 MCMC 方法，但我从未找到任何最佳值，主要是因为问题 2

2) 步骤 2-4 本身就非常低效。

我尝试过类似的东西（类似于伪代码；实际的功能是组成的）。关于在 A、B 上使用广播可能存在许多小问题，但这些问题不如需要为每一步进行插值的问题重要。

我认识的人推荐使用某种Expectation Maximization 算法，但我对此知之甚少，无法从头开始编写代码。我真的希望有一些很棒的 scipy（或其他开源）算法我无法找到解决我的整个问题，但在这一点上我不抱希望。

import numpy as np
import scipy as sci
from scipy import interpolate

X_data
Y_data

def x(t,A,B):
    return A**t + B**t
def y(t,A,B):
    return A*t + B

def interp(A,B):
    ts = np.arange(-10,10,0.1)
    xs = x(ts,A,B)
    ys = y(ts,A,B)
    f = interpolate.interp1d(xs,ys)
    return f

N = 101
lsqs = np.recarray((N**2),dtype=float)

count = 0
for i in range(0,N):
    A = 0.1*i            #checks A between 0 and 10
    for j in range(0,N):
        B = 10 + 0.1*j   #checks B between 10 and 20

        f = interp(A,B)
        y_fit = f(X_data)
        squares = np.sum((y_fit - Y_data)**2)

        lsqs[count] = (A,b,squares) #puts the values in place for comparison later
        count += 1        #allows us to move to the next cell

i = np.argmin(lsqs[:,2])

A_optimal = lsqs[i][0]
B_optimal = lsqs[i][1]

【问题讨论】：

你应该看看scipy.interpolate.splprep

不幸的是，这仅对连接我的(X,Y) 数据上的点有用，不适用于查找*params

【参考方案1】：

如果我正确理解了这个问题，参数是每个样本中相同的常量，但 t 因样本而异。因此，例如，也许您有一大堆点，您认为这些点是从一个圆圈中采样的

x = a+r cos(t)   
y = b+r sin(t)

t 的不同值。

在这种情况下，我要做的是消除变量t 以获得x 和y 之间的关系——在这种情况下，(x-a)^2+(y-b)^2 = r^2。如果您的数据完全符合模型，您的每个数据点都会有(x-a)^2+(y-b)^2 = r^2。有一些错误，你仍然可以找到(a,b,r) 来最小化

sum_i ((x_i-a)^2 + (y_i-b)^2 - r^2)^2.

Mathematica 的Eliminate 命令在某些情况下可以自动执行消除 t 的过程。

PS 你可能在 stats.stackexchange、math.stackexchange 或 mathoverflow.net 上做得更好。我知道最后一个的名声很吓人，但我们不咬人，真的！

【讨论】：

感谢您的反馈。我正在尝试完全按照您的建议进行操作 - 找到 a 和 b 的等价物。然而，方程并不总是解析求解 t，当它们解析求解时，它们最终会产生许多奇点。当我按照t 保留它时，我只有t=0 的奇点。但是，当我消除 t 时，我得到的部分方程看起来像（为简单起见，去掉了一些常数）y - x + c' = y/x + (y/(x-y))^(1/3) 这意味着 x 中相对较小的错误可能会压倒拟合过程。