使用Scipy Optimize库拟合数据时不稳定

klr1opcd  于 2022-11-10  发布在  其他
关注(0)|答案(1)|浏览(186)

我正在尝试用一个函数来拟合我所拥有的一组数据。这个函数是:
x(t)= - B +平方(AB(t-t0)+(x 0 + B)^2)
我试着把我的数据(包含在底部),但我发现无论我做什么,B的拟合都是极不稳定的。改变方法或初始猜测都会极大地改变输出值。此外,当我使用curve_fit查看此拟合的误差时,误差几乎比值高两个数量级。有人对我应该做些什么来减少错误有什么建议吗?

import numpy as np
import scipy.optimize as spo

def modelFun(t,A,B):
    return -B + np.sqrt(A*B*(t-t0) + np.power(x0 + B,2))

def errorFun(k,time,data):
    A = k[0]
    B = k[1]
    return np.sum((data-modelFun(time,A,B))**2)

data = np.genfromtxt('testdata.csv',delimiter=',',skip_header = 1)
time = data[:,0]
xt = data[:,1]
t0 = data[0,0]
x0 = data[0,1]

minErrOut = spo.minimize(errorFun,[1,1000],args=(time,xt),bounds=((0,None),(0,None)))
(curveOut, curveCovar) = spo.curve_fit(modelFun,time,xt,p0=[1,1000],method='dogbox',bounds=([-np.inf,0],np.inf))

print('minimize result: A={}; B={}'.format(*minErrOut.x))
print('curveFit result: A={}; B={}'.format(*curveOut))
print('curveFit Error: A={}; B={}'.format(*np.sqrt(np.diag(curveCovar))))

数据文件:

Time,x
201,2.67662
204,3.28159
206,3.44378
208,3.72537
210,3.94826
212,4.36716
214,4.65373
216,5.26766
219,5.59502
221,6
223,6.22189
225,6.49652
227,6.799
229,7.30846
231,7.54229
232,7.76517
233,7.6209
234,7.89552
235,7.94826
236,8.17015
237,8.66965
238,8.66965
239,8.8398
240,8.88856
241,9.00697
242,9.45075
243,9.51642
244,9.63483
245,9.63483
246,10.07861
247,10.02687
248,10.24876
249,10.31443
250,10.47164
251,10.99502
252,10.92935
253,11.0995
254,11.28358
255,11.58209
256,11.53035
257,11.62388
258,11.93632
259,11.98806
260,12.26269
261,12.43284
262,12.60299
263,12.801
264,12.99502
265,13.08557
266,13.25572
267,13.32139
268,13.57114
269,13.76617
270,13.88358
271,13.83184
272,14.10647
273,14.27662
274,14.40796
iecba09b

iecba09b1#

TL;DR;

您的数据集是线性的,在较大的时间尺度上会丢失观测值。因此,您可以捕获与斜率成比例的A,而您的模型需要使B保持较大(可能是无界的)以抑制平方根趋势。
这可以通过展开模型的泰勒级数并分析与回归关联的MSE曲面来确认。
简单地说,考虑到这类数据集和给定的模型,接受A不信任B
MCVE公司
首先,让我们重现您的问题:

import io
import numpy as np
import pandas as pd
from scipy import optimize

stream = io.StringIO("""Time,x
201,2.67662
204,3.28159
...
273,14.27662
274,14.40796""")
data = pd.read_csv(stream)

# Origin Shift:

data = data.sub(data.iloc[0,:])
data = data.set_index("Time")

# Simplified model:

def model(t, A, B):
    return -B + np.sqrt(A*B*t + np.power(B, 2))

# NLLS Fit:

parameters, covariance = optimize.curve_fit(model, data.index, data["x"].values, p0=(1, 1000), ftol=1e-8)

# array([3.23405915e-01, 1.59960168e+05])

# array([[ 3.65068730e-07, -3.93410484e+01],

# [-3.93410484e+01,  9.77198860e+12]])

调整是公平的:

但是,正如您所注意到的,模型参数相差几个数量级,这可能会妨碍优化正常执行。
请注意,您的数据集是非常线性的。观察到的效应并不令人惊讶,并且是所选模型固有的。B参数必须比A大几个数量级才能保持线性行为。
first terms of the Taylor series的分析支持这一说法:

def taylor(t, A, B):
    return (A/2*t - A**2/B*t**2/8) 

parameters, covariance = optimize.curve_fit(taylor, data.index, data["x"].values, p0=(1, 100), ftol=1e-8)
parameters

# array([3.23396685e-01, 1.05237134e+09])

毫无疑问,可以捕获线性数据集的斜率,而参数B可以任意大,并将在优化过程中导致浮点运算错误(因此您会看到下面的minimize警告)。

分析误差曲面

这个问题可以重新表述为一个最小化问题:

def objective(beta, t, x):
    return np.sum(np.power(model(t, beta[0], beta[1]) - x, 2))

result = optimize.minimize(objective, (1, 100), args=(data.index, data["x"].values))

# fun: 0.6594398116927569

# hess_inv: array([[8.03062155e-06, 2.94644208e-04],

# [2.94644208e-04, 1.14979735e-02]])

# jac: array([2.07304955e-04, 6.40749931e-07])

# message: 'Desired error not necessarily achieved due to precision loss.'

# nfev: 389

# nit: 50

# njev: 126

# status: 2

# success: False

# x: array([3.24090627e-01, 2.11891188e+03])

如果我们绘制与数据集关联的MSE,则会得到以下曲面:

我们在A空间上有一个狭窄的峡谷,但在B空间上至少在最初的几十年里似乎是无界的。这支持了你在帖子和评论中的观察。它也带来了一个技术上的见解,为什么我们不能正确地适应B
对合成数据集执行相同的操作:

t = np.linspace(0, 1000, 100)
x = model(t, 0.35, 20)
data = pd.DataFrame(x, index=t, columns=["x"])

除了开始时的线性趋势外,还具有平方根形状。

result = optimize.minimize(objective, (1, 0), args=(data.index, data["x"].values), tol=1e-8)

# fun: 1.9284246829733202e-10

# hess_inv: array([[ 4.34760333e-05, -4.42855253e-03],

# [-4.42855253e-03,  4.59219063e-01]])

# jac: array([ 4.35726463e-03, -2.19158602e-05])

# message: 'Desired error not necessarily achieved due to precision loss.'

# nfev: 402

# nit: 94

# njev: 130

# status: 2

# success: False

# x: array([ 0.34999987, 20.000013  ])

此版本的问题具有以下MSE表面:
第一次
在已知解周围显示潜在的凸谷,这解释了当有足够长的时间采集时,为什么你能够拟合两个参数。
请注意,谷值被强烈拉伸,这意味着在这种情况下,您的问题将受益于规范化。

相关问题