我有 (X,Y)
形式的实验数据和 (x(t;*params),y(t;*params))
其中 t
是物理(但不可观察)变量,*params
是我要确定的参数。 t
是一个连续变量,x
和t
、y
和y
之间是1:1的关系t
在模型中。
在一个完美的世界中,我会知道 T
的值(参数的实际值),并且能够进行极其基本的最小二乘法拟合以找到*参数
。 (请注意,我不试图在我的情节中“连接”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]
最佳答案
如果我正确理解了这个问题,参数是每个样本中都相同的常量,但 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 上做得更好。我知道最后一个名声不好,但我们不咬人,真的!
关于python - 在 Python 中拟合参数曲线,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/32133733/