终身教育平台
5.2.3 多项式回归分析
理论上可以利用上小节的方法拟合任何曲线，但前提条件是要事先对模型的形式进行判断，即知道非线性模型的参数设置。在一般情况下，通过绘制散点图的形式可以做到这一点
但是在更一般的情况下，如有多个自变量的情况下，无法绘制散点图，同时也很难对模型形式进行预估，这个时候可以使用本小节所介绍的方法。根据数学的相关理论，任何曲线均可以使用多项式进行逼近，这种逼近的分析过程即多项式回归。
多项式回归类似于可线性化的非线性模型，可通过变量代换的方式使用普通最小二乘对参数进行估计。
设有因变量y和自变量x，它们之间的关系为n次多项式的关系，则有如下模型：
令x1=x,x2=x2,…,xn=xn，则多项式模型就转化为如下的多元线性模型：
这样就可以按照多元线性回归模型进行分析了。对于多元的多项式模型，同样可以进行变量代换，转化之后的模型同样可以按照多元线性回归模型进行分析。多项式回归当阶数过高时，待估参数过多，在样本量不大的情况下会比较困难，这是多项式回归的一大缺陷。因此，一般的多项式回归模型很少应用到三阶以上。
例5-4
基尼系数是20世纪初意大利学者基尼根据洛仑兹曲线所定义用于判断收入分配公平程度的指标。其值在0和1之间，越往1靠近则表示收入分配公平程度越低。如果要较为精确的计算基尼系数，则主要通过拟合洛仑兹曲线来进行。其中拟合洛仑兹曲线的方法主要有广义二次洛仑兹曲线法、洛伦兹曲线法等，本书利用CHNS（中国营养调查）中2008年的人均家庭收入数据，按照1322个有效调查样本分别计算出了收入和人口的累计百分比，如“lorenz.csv”所示。试利用本例数据对洛仑兹曲线进行拟合。
In [34]:
lorenz=pd.read_csv('lorenz.csv')
lorenz.head()
Out [34]:
| cpop | cincome | |
| 0 | 0.018 | 0.000000e+00 |
| 1 | 0.019 | 9.701260e-07 |
| 2 | 0.020 | 3.366135e-06 |
| 3 | 0.020 | 5.785387e-06 |
| 4 | 0.021 | 8.521446e-06 |
由于事先不知道曲线的形式，因此无法像非线性回归模型那样为模型设置参数，故考虑采用多项式回归分析。numpy.polyfit可以直接对多项式函数进行拟合：
In [35]:
lorenz_est=np.polyfit(lorenz['cpop'],lorenz['cincome'],2)
#本例中三个参数分别是自变量、因变量和阶数
lorenz_est
Out[35]:
array([ 1.02722918, -0.20824691, 0.02550764])
估计出参数之后，可以使用numpy的poly1d对象来处理多项式，此对象可以像函数一样调用，并且返回多项式的值：
In [36]:
lorenz_poly=np.poly1d(lorenz_est)
print lorenz_poly
Out [36]:
2
1.027 x - 0.2082 x + 0.02551
上述结果即为估计出来的多项式回归方程式。
在实际应用中，除了得到上述拟合的回归方程式之外，更多的是想要得到模型的检验过程。如前所述，多项式回归模型可通过转换之后按照多元线性回归模型进行分析。因此，同样可以使用statsmodels中的ols方法进行分析：
In [37]:
formula='cincome~cpop+np.square(cpop)'
lorenz_model1=ols(formula,data=lorenz).fit()
print lorenz_model1.summary2()
Out [37]:
Results: Ordinary least squares
===================================================================
Model: OLS Adj. R-squared: 0.994
Dependent Variable: cincome AIC: -6659.2479
Date: 2017-09-09 16:08 BIC: -6643.6872
Df Model: 2 F-statistic: 1.105e+05
Df Residuals: 1319 Prob (F-statistic): 0.00
R-squared: 0.994 Scale: 0.00037922
R-squared: 0.994 Scale: 0.00037922
-------------------------------------------------------------------
Coef. Std.Err. t P>|t| [0.025 0.975]
-------------------------------------------------------------------
Intercept 0.0255 0.0017 14.9508 0.0000 0.0222 0.0289
cpop -0.2082 0.0077 -26.9092 0.0000 -0.2234 -0.1931
np.square(cpop) 1.0272 0.0074 139.4493 0.0000 1.0128 1.0417
-------------------------------------------------------------------
Omnibus: 732.937 Durbin-Watson: 0.002
Prob(Omnibus): 0.000 Jarque-Bera (JB): 8902.124
Skew: 2.320 Prob(JB): 0.000
Kurtosis: 14.835 Condition No.: 24
===================================================================
在上述结果可以看到，cpop和square(cpop)即二次项在模型中均显著，且F统计量为1.105e+05，模型非常显著。模型的R2也非常高，达到了0.9941且是非常显著的。此外，参数估计和检验结果表明，各项回归系数均非常显著。据此可写出形如lorenz_poly对象的多项式回归方程（本书不再赘述）。
In [38]:
plot_fit(lorenz_model1,1)
plt.savefig('9278.svg',bbox_inches='tight')
plt.show()
Out [38]:
In [39]:
plt.plot(lorenz['cpop'],lorenz['cincome'],'ro',lorenz['cpop'],lorenz_model1.fittedvalues,linewidth=2)
Out[39]:
[<matplotlib.lines.Line2D at 0x1135dd950>,
<matplotlib.lines.Line2D at 0x1135dd050>]
由上图可以看出，本例所拟合出来的二次曲线与原始数据散点图较为接近，只是在人口累计百分比的两侧，偏离程度有点大。
但对于从洛仑兹曲线中进行基尼系数测算的这个实际问题而言，要在准确拟合曲线的基础上才能精确的计算出基尼系数。基于此，本例可以考虑使用高次（三次）多项式重新拟合曲线：
In [40]:
#因为statsmodels的ols方法建立模型时变量不能使用表达式做变换，只能采用函数的方式进行变换，故我们先定义一个三次方函数：
def cube(x):
return x**3
formula='cincome~cpop+np.square(cpop)+cube(cpop)'
lorenz_model2=ols(formula,data=lorenz).fit()
print lorenz_model2.summary2()
Out [40]:
Results: Ordinary least squares
===================================================================
Model: OLS Adj. R-squared: 0.998
Dependent Variable: cincome AIC: -7995.1774
Date: 2017-07-20 16:47 BIC: -7974.4298
No. Observations: 1322 Log-Likelihood: 4001.6
Df Model: 3 F-statistic: 2.034e+05
Df Residuals: 1318 Prob (F-statistic): 0.00
R-squared: 0.998 Scale: 0.00013794
-------------------------------------------------------------------
Coef. Std.Err. t P>|t| [0.025 0.975]
-------------------------------------------------------------------
Intercept -0.0238 0.0015 -16.4017 0.0000 -0.0267 -0.0210
cpop 0.3350 0.0122 27.3831 0.0000 0.3110 0.3589
np.square(cpop) -0.2877 0.0277 -10.3769 0.0000 -0.3421 -0.2333
cube(cpop) 0.8611 0.0179 48.0430 0.0000 0.8259 0.8963
-------------------------------------------------------------------
Omnibus: 943.668 Durbin-Watson: 0.005
Prob(Omnibus): 0.000 Jarque-Bera (JB): 23537.954
Skew: 2.997 Prob(JB): 0.000
Kurtosis: 22.784 Condition No.: 134
===================================================================
从上述结果中可以看到，回归方程拟合程度的R2更高，回归方程总体显著性检验的F值为2.034e+05，对应P值几乎为0，非常显著。参数估计及回归系数显著性检验结果表明，各回归系数均显著。因此，可以写出如下的高次多项式回归方程：
cincome=-0.0238+0.3350×cpop-0.2877×cpop2+0.8611×cpop3
最后修改: 2020年07月24日 Friday 09:40
终身教育平台