在Scipy中曲线拟合真实数据时如何修复“RuntimeWarning：exp遇到溢出”？

Question

我正在尝试确定用于现实世界数据的最佳模型函数和参数。我有几个数据集，它们都表现出类似的指数衰减，我想计算每个数据集的拟合函数参数。

最大的数据集在 x 轴上从 1 到大约 1,000,000 不等，在 y 轴上从 0 到大约 10,000 不等。

我是 Numpy 和 Scipy 的新手，所以我尝试将这个问题的代码适应我的数据，但没有成功： 在没有初始猜测的情况下拟合指数衰减

# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
import scipy as sp
import scipy.optimize

x = np.array([   1.,    4.,    9.,   16.,   25.,   36.,   49.,   64.,   81.,  100.,  121.,
        144.,  169.,  196.,  225.,  256.,  289.,  324.,  361.,  400.,  441.,  484.,
        529.,  576.,  625.,  676.,  729.,  784.,  841.,  900.,  961., 1024., 1089.,
       1156., 1225., 1296., 1369., 1444., 1521., 1600., 1681., 1764., 1849., 1936.,
       2025., 2116., 2209., 2304., 2401., 2500., 2601., 2704., 2809., 2916., 3025.,
       3136., 3249., 3364., 3481., 3600., 3721., 3844., 3969., 4096., 4225., 4356.,
       4489., 4624., 4761., 4943.])

y = np.array([3630., 2590., 2063., 1726., 1484., 1301., 1155., 1036.,  936.,  851.,  778.,
        714.,  657.,  607.,  562.,  521.,  485.,  451.,  421.,  390.,  362.,  336.,
        312.,  293.,  279.,  265.,  253.,  241.,  230.,  219.,  209.,  195.,  183.,
        171.,  160.,  150.,  142.,  134.,  127.,  120.,  114.,  108.,  102.,   97.,
         91.,   87.,   83.,   80.,   76.,   73.,   70.,   67.,   64.,   61.,   59.,
         56.,   54.,   51.,   49.,   47.,   45.,   43.,   41.,   40.,   38.,   36.,
         35.,   33.,   31.,   30.])

# Define the model function
def model_func(x, A, K, C):
    return A * np.exp(-K * x) + C

# Optimise the curve
opt_parms, parm_cov = sp.optimize.curve_fit(model_func, x, y, maxfev=1000)

# Fit the parameters to the data
A, K, C = opt_parms
fit_y = model_func(x, A, K, C)

# Visualise the original data and the fitted function
plt.clf()
plt.title('Decay Data')
plt.plot(x, y, 'r.', label='Actual Data\n')
plt.plot(x, fit_y, 'b-', label='Fitted Function:\n $y = %0.2f e^{%0.2f t} + %0.2f$' % (A, K, C))
plt.legend(bbox_to_anchor=(1, 1), fancybox=True, shadow=True)
plt.show()

当我使用 Python 2.7（在 Windows 7 64 位上）运行此代码时，我收到错误消息RuntimeWarning: overflow encountered in exp。上图显示了该函数不适合我的数据的问题。

我使用的模型函数对我的数据正确吗？如果是这样，我怎样才能更好地计算拟合函数参数，以便我可以将它们与新数据一起使用？

Answer 1

我认为你需要记录 X 数据。这是我的二次对数方程“y = a + b * np.log(x) + c * np.log(x) ** 2”的图。所有 1.0 的 scipy 默认初始参数估计对我来说效果很好，给出 RMSE = 4.020 和 R 平方 = 0.9999。

更新：添加了 Python 图形拟合器

import numpy, scipy, matplotlib
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

xData = numpy.array([   1,    4,    9,   16,   25,   36,   49,   64,   81,  100,  121,
        144,  169,  196,  225,  256,  289,  324,  361,  400,  441,  484,
        529,  576,  625,  676,  729,  784,  841,  900,  961, 1024, 1089,
       1156, 1225, 1296, 1369, 1444, 1521, 1600, 1681, 1764, 1849, 1936,
       2025, 2116, 2209, 2304, 2401, 2500, 2601, 2704, 2809, 2916, 3025,
       3136, 3249, 3364, 3481, 3600, 3721, 3844, 3969, 4096, 4225, 4356,
       4489, 4624, 4761, 4943], dtype=float)

yData = numpy.array([3630, 2590, 2063, 1726, 1484, 1301, 1155, 1036,  936,  851,  778,
        714,  657,  607,  562,  521,  485,  451,  421,  390,  362,  336,
        312,  293,  279,  265,  253,  241,  230,  219,  209,  195,  183,
        171,  160,  150,  142,  134,  127,  120,  114,  108,  102,   97,
         91,   87,   83,   80,   76,   73,   70,   67,   64,   61,   59,
         56,   54,   51,   49,   47,   45,   43,   41,   40,   38,   36,
         35,   33,   31,   30], dtype=float)


def func(x, a, b, c): # simple quadratic example
    return a + b*numpy.log(x) + c*numpy.log(x)**2


# these are the same as the scipy defaults
initialParameters = numpy.array([1.0, 1.0, 1.0])

# curve fit the test data
fittedParameters, pcov = curve_fit(func, xData, yData, initialParameters)

modelPredictions = func(xData, *fittedParameters) 

absError = modelPredictions - yData

SE = numpy.square(absError) # squared errors
MSE = numpy.mean(SE) # mean squared errors
RMSE = numpy.sqrt(MSE) # Root Mean Squared Error, RMSE
Rsquared = 1.0 - (numpy.var(absError) / numpy.var(yData))

print('Parameters:', fittedParameters)
print('RMSE:', RMSE)
print('R-squared:', Rsquared)

print()


##########################################################
# graphics output section
def ModelAndScatterPlot(graphWidth, graphHeight):
    f = plt.figure(figsize=(graphWidth/100.0, graphHeight/100.0), dpi=100)
    axes = f.add_subplot(111)

    # first the raw data as a scatter plot
    axes.plot(xData, yData,  'D')

    # create data for the fitted equation plot
    xModel = numpy.linspace(min(xData), max(xData))
    yModel = func(xModel, *fittedParameters)

    # now the model as a line plot
    axes.plot(xModel, yModel)

    axes.set_xlabel('X Data') # X axis data label
    axes.set_ylabel('Y Data') # Y axis data label

    plt.show()
    plt.close('all') # clean up after using pyplot

graphWidth = 800
graphHeight = 600
ModelAndScatterPlot(graphWidth, graphHeight)