我有形式的实验数据和形式(X,Y)
的理论模型,(x(t;*params),y(t;*params))
其中t
是物理(但不可观察)变量,并且*params
是我想要确定的参数。是一个连续变量,在模型中和和之间t
存在1:1的关系。x
t
y
t
在一个完美的世界里,我会知道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 上使用广播可以提出许多小问题,但这些问题不如需要为每一步进行插值的问题重要。
我认识的人建议使用某种期望最大化算法,但我对此知之甚少,无法从头开始编写代码。我真的希望有一些很棒的 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]