为什么 statsmodels 的 OLS 中的四次线性回归与 LibreOffice Calc 不匹配？答案

【问题标题】：Why doesn't a quartic linear regression in statsmodels' OLS match LibreOffice Calc?为什么 statsmodels 的 OLS 中的四次线性回归与 LibreOffice Calc 不匹配？
【发布时间】：2019-07-03 01:31:12
【问题描述】：

我正在使用 statsmodels 的 OLS 线性回归和 Patsy 四次公式 y ~ x + I(x**2) + I(x**3) + I(x**4)，但与 LibreOffice Calc 相比，所得回归与数据的拟合度较差。为什么这与 LibreOffice Calc 生成的不匹配？

statsmodels 代码：

import io
import numpy
import pandas
import matplotlib
import matplotlib.offsetbox
import statsmodels.tools
import statsmodels.formula.api

csv_data = """Year,CrudeRate
1999,197.0
2000,196.5
2001,194.3
2002,193.7
2003,192.0
2004,189.2
2005,189.3
2006,187.6
2007,186.9
2008,186.0
2009,185.0
2010,186.2
2011,185.1
2012,185.6
2013,185.0
2014,185.6
2015,185.4
2016,185.1
2017,183.9
"""

df = pandas.read_csv(io.StringIO(csv_data))

cause = "Malignant neoplasms"
x = df["Year"].values
y = df["CrudeRate"].values

olsdata = {"x": x, "y": y}
formula = "y ~ x + I(x**2) + I(x**3) + I(x**4)"
model = statsmodels.formula.api.ols(formula, olsdata).fit()

print(model.params)

df.plot("Year", "CrudeRate", kind="scatter", grid=True, title="Deaths from {}".format(cause))

func = numpy.poly1d(model.params.values[::-1])
matplotlib.pyplot.plot(df["Year"], func(df["Year"]))

matplotlib.pyplot.show()

产生以下系数：

Intercept    9.091650e-08
x            9.127904e-05
I(x ** 2)    6.109623e-02
I(x ** 3)   -6.059164e-05
I(x ** 4)    1.503399e-08

还有下图：

但是，如果我将数据带入 LibreOffice Calc，单击绘图并选择“插入趋势线...”，选择“多项式”，输入“度数”=4，然后选择“显示方程”，结果趋势线与 statsmodels 不同，似乎更适合：

系数是：

Intercept = 1.35e10
x =          2.69e7
x^2 =       -2.01e4
x^3 =          6.69
x^4 =      -0.83e-3

statsmodels 版本：

$ pip3 list | grep statsmodels
statsmodels                  0.9.0

编辑：三次也不匹配，但二次匹配。

编辑：缩小Year（并在 LibreOffice 中做同样的事情）匹配：

df = pandas.read_csv(io.StringIO(csv_data))
df["Year"] = df["Year"] - 1998

按比例缩小后的系数和绘图：

Intercept    197.762384
x             -0.311548
I(x ** 2)     -0.315944
I(x ** 3)      0.031304
I(x ** 4)     -0.000833

【问题讨论】：

我的猜测是回归中的 X 矩阵的条件很差，因为年份值很大。试试year - 1998 作为趋势变量。
也可以按比例缩小，x**4 相对于常数 1 来说会非常大。
@JamesPhillips 三阶立方也与 LibreOffice 完全不同。如果我下降到二阶二次，那么事情就会匹配。
statsmodels 不会进行任何自动重新缩放。多项式不适用于大数，应始终缩放到“合理”范围。例如 numpy.polynomial 可以选择缩放到所有多项式都表现良好的区间 [-1, 1]。
@Josef 为了澄清我的最后一条评论，如果我通过减去 1998 在 LibreOffice 中执行相同的缩放，那么所有系数都匹配并且我的问题得到解决，所以我只剩下问题了关于我是否应该总是缩小到个位数？

标签： python pandas curve-fitting statsmodels

【解决方案1】：

基于来自@Josef 的 cmets，问题在于大数不适用于高阶多项式，并且 statsmodels 无法自动缩放域。另外，我在原始问题中没有提到这一点，因为我没想到需要转换域，但我还需要根据年份预测样本外值，所以我将其设为范围结束：

predict_x = +5
min_scaled_domain = -1
max_scaled_domain = +1
df["Year"] = df["Year"].transform(lambda x: numpy.interp(x, (x.min(), x.max() + predict_x), (min_scaled_domain, max_scaled_domain)))

这种转换创建了一个拟合良好的回归：

如果在 LibreOffice Calc 中应用相同的域变换，则系数匹配。

最后，打印预测值：

func = numpy.polynomial.Polynomial(model.params)
print(func(max_scaled_domain))

【讨论】：

我将未转换的原始数据很好地拟合为带偏移量的标准几何方程“y = a * pow(x, (b * x)) + Offset”，它有两个形状参数“a”和“b”加上不影响曲线形状的偏移项。拟合参数 a = 1.0774214144669925E+131，b = -1.9682480254878586E-02，Offset = 1.8369038653939549E+02 得出 R-squared = 0.968 和 RMSE = 0.734
@JamesPhillips 我不关注 - 你是如何确定 a、b 和 Offset 的？
既然你问了，我在我的开源 zunzun.com 在线曲线拟合网站上使用了“函数查找器”来将数据拟合到一百多个已知的命名方程，这个看起来像一个极好的候选人。然后，我将您发布的数据粘贴到 zunzun.com/Equation/2/Power/… 的该方程式的“文本数据编辑器”中，然后点击提交按钮。