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

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

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

  1. import numpy as np
  2. import scipy.optimize as spo
  3. def modelFun(t,A,B):
  4. return -B + np.sqrt(A*B*(t-t0) + np.power(x0 + B,2))
  5. def errorFun(k,time,data):
  6. A = k[0]
  7. B = k[1]
  8. return np.sum((data-modelFun(time,A,B))**2)
  9. data = np.genfromtxt('testdata.csv',delimiter=',',skip_header = 1)
  10. time = data[:,0]
  11. xt = data[:,1]
  12. t0 = data[0,0]
  13. x0 = data[0,1]
  14. minErrOut = spo.minimize(errorFun,[1,1000],args=(time,xt),bounds=((0,None),(0,None)))
  15. (curveOut, curveCovar) = spo.curve_fit(modelFun,time,xt,p0=[1,1000],method='dogbox',bounds=([-np.inf,0],np.inf))
  16. print('minimize result: A={}; B={}'.format(*minErrOut.x))
  17. print('curveFit result: A={}; B={}'.format(*curveOut))
  18. print('curveFit Error: A={}; B={}'.format(*np.sqrt(np.diag(curveCovar))))

数据文件:

  1. Time,x
  2. 201,2.67662
  3. 204,3.28159
  4. 206,3.44378
  5. 208,3.72537
  6. 210,3.94826
  7. 212,4.36716
  8. 214,4.65373
  9. 216,5.26766
  10. 219,5.59502
  11. 221,6
  12. 223,6.22189
  13. 225,6.49652
  14. 227,6.799
  15. 229,7.30846
  16. 231,7.54229
  17. 232,7.76517
  18. 233,7.6209
  19. 234,7.89552
  20. 235,7.94826
  21. 236,8.17015
  22. 237,8.66965
  23. 238,8.66965
  24. 239,8.8398
  25. 240,8.88856
  26. 241,9.00697
  27. 242,9.45075
  28. 243,9.51642
  29. 244,9.63483
  30. 245,9.63483
  31. 246,10.07861
  32. 247,10.02687
  33. 248,10.24876
  34. 249,10.31443
  35. 250,10.47164
  36. 251,10.99502
  37. 252,10.92935
  38. 253,11.0995
  39. 254,11.28358
  40. 255,11.58209
  41. 256,11.53035
  42. 257,11.62388
  43. 258,11.93632
  44. 259,11.98806
  45. 260,12.26269
  46. 261,12.43284
  47. 262,12.60299
  48. 263,12.801
  49. 264,12.99502
  50. 265,13.08557
  51. 266,13.25572
  52. 267,13.32139
  53. 268,13.57114
  54. 269,13.76617
  55. 270,13.88358
  56. 271,13.83184
  57. 272,14.10647
  58. 273,14.27662
  59. 274,14.40796
iecba09b

iecba09b1#

TL;DR;

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

  1. import io
  2. import numpy as np
  3. import pandas as pd
  4. from scipy import optimize
  5. stream = io.StringIO("""Time,x
  6. 201,2.67662
  7. 204,3.28159
  8. ...
  9. 273,14.27662
  10. 274,14.40796""")
  11. data = pd.read_csv(stream)
  12. # Origin Shift:
  13. data = data.sub(data.iloc[0,:])
  14. data = data.set_index("Time")
  15. # Simplified model:
  16. def model(t, A, B):
  17. return -B + np.sqrt(A*B*t + np.power(B, 2))
  18. # NLLS Fit:
  19. parameters, covariance = optimize.curve_fit(model, data.index, data["x"].values, p0=(1, 1000), ftol=1e-8)
  20. # array([3.23405915e-01, 1.59960168e+05])
  21. # array([[ 3.65068730e-07, -3.93410484e+01],
  22. # [-3.93410484e+01, 9.77198860e+12]])

调整是公平的:

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

  1. def taylor(t, A, B):
  2. return (A/2*t - A**2/B*t**2/8)
  3. parameters, covariance = optimize.curve_fit(taylor, data.index, data["x"].values, p0=(1, 100), ftol=1e-8)
  4. parameters
  5. # array([3.23396685e-01, 1.05237134e+09])

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

分析误差曲面

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

  1. def objective(beta, t, x):
  2. return np.sum(np.power(model(t, beta[0], beta[1]) - x, 2))
  3. result = optimize.minimize(objective, (1, 100), args=(data.index, data["x"].values))
  4. # fun: 0.6594398116927569
  5. # hess_inv: array([[8.03062155e-06, 2.94644208e-04],
  6. # [2.94644208e-04, 1.14979735e-02]])
  7. # jac: array([2.07304955e-04, 6.40749931e-07])
  8. # message: 'Desired error not necessarily achieved due to precision loss.'
  9. # nfev: 389
  10. # nit: 50
  11. # njev: 126
  12. # status: 2
  13. # success: False
  14. # x: array([3.24090627e-01, 2.11891188e+03])

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

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

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

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

  1. result = optimize.minimize(objective, (1, 0), args=(data.index, data["x"].values), tol=1e-8)
  2. # fun: 1.9284246829733202e-10
  3. # hess_inv: array([[ 4.34760333e-05, -4.42855253e-03],
  4. # [-4.42855253e-03, 4.59219063e-01]])
  5. # jac: array([ 4.35726463e-03, -2.19158602e-05])
  6. # message: 'Desired error not necessarily achieved due to precision loss.'
  7. # nfev: 402
  8. # nit: 94
  9. # njev: 130
  10. # status: 2
  11. # success: False
  12. # x: array([ 0.34999987, 20.000013 ])

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

展开查看全部

相关问题