在 Python 中拟合参数曲线

Nec*_*ion 7 python model-fitting parametric-equations

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

在完美的世界中,我会知道(参数的真实世界值)的值T,并且能够进行极其基本的最小二乘拟合来找到 的值*params。(请注意,我并没有x尝试在图中“连接”和的值y,例如3124300231464345。)我无法保证在我的实际数据中,潜在值T是单调的,因为我的数据是跨多个周期收集的。

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

  1. 选择 的一些值*params并将其应用于模型
  2. 获取一个值数组t并将其放入模型中以创建一个数组model(*params) = (x(*params),y(*params))
  3. X(数据值)插值model得到Y_predicted
  4. Y在和之间运行最小二乘(或其他)比较Y_predicted
  5. 再次执行此操作以获得一组新的*params
  6. 最终选择最佳值*params

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

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

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

我尝试过类似的东西(类似于伪代码;实际的函数是组成的)。关于在 A、B 上使用广播可能会出现许多小问题,但这些问题没有需要对每个步骤进行插值的问题那么重要。

我认识的人建议使用某种期望最大化算法,但我对此了解不够,无法从头开始编写代码。我真的希望有一些很棒的 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]
Run Code Online (Sandbox Code Playgroud)

小智 0

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

x = a+r cos(t)   
y = b+r sin(t)
Run Code Online (Sandbox Code Playgroud)

在不同的值t.

在这种情况下,我要做的就是消除变量以获得和t之间的关系——在这种情况下,是 。如果您的数据完全符合模型,那么您的每个数据点都会有。如果出现一些错误,您仍然可以发现最小化xy(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 上做得更好。我知道最后一种有可怕的名声,但我们不咬人,真的!