利用python进行回归分析

通常大家会认为曲线拟合和回归分析类似，但其实回归分析中是包含曲线拟合的。拟合是研究因变量和自变量的函数关系的。而回归是研究随机变量间的相关关系的。拟合侧重于调整参数，使得与给出的数据相符合。而回归则是侧重于研究变量的关系，对拟合问题做统计分析。

一元线性回归

一元线性回归模型的一般形式

数据通常呈一条直线，则y和x之间的关系通常可以看做近似线性关系。但是一般来说这些数据点并不在一条直线上，这说明y和x的关系并没有确切到给定x就可以唯一确定y的程度。其实y还受到很多因素的影响。如果主要研究y和x的关系，可以假设有如下关系。

（1）为未知待定常数称为回归系数，是其他随机因素对y的影响，并且服从分布。

（2）对于参数的估计有多重方法，一方面我们可以使用拟合的方法去估计参数，另一方面我们可以用最小二乘法进行参数的估计。

最小二乘法：使得回归模型和直线方程所在的所有数据点都比较接近。

,我们为了求得最小Q，对这个式子求导最终计算出。

，，，

（3）相关性检验，拟合优度和剩余标准差:

,这是原始数据的总变异平方和。

,拟合直线可解释的变异平方和。

,残差平方和。

由公式推导可以证明有：

对于一个具体的样本来说，是一个定值。所以可解释变异越大，则必有残差平方和越小。这个式子可以从两个方面解释拟合程度的优良。

（1）越大，在回归方程中解释因变量得变异部分越大，回归方程对原始数据的解释越好。

（2）越小，观测值围绕回归直线越紧密，回归方程对源数据的拟合程度越好，

 因此定义一个测量标量来说明回归方程对原始数据的拟合程度，这个系数叫拟合优度。拟合优度是指可解释变异占总体变异的百分比，用表示。有：

有以下简单的性质：

（1）

（2）当=1时，有,此时原数据的变异程度完全可以由拟合值的变异来解释，且残差为0（）,即原数据和拟合点完全吻合。

（3） 当=0时，此时的回归方程完全不能解释原数据的总变异，且残差为0（）,即y得变异和x完全无关。

 （ 4）回归方程的显著检验

在以上的讨论中，我们假设y是关于x的回归方程具有形式。在实际中，需要检验是否为x的线性函数，若与x成线性函数为真，则 不为0，如果为0，则y和x就无线性关系，所以我们需要进行假设检验。 

提出假设：:,：

如果接受原假设，则说明回归结果不显著，否则回归结果显著。

求解可到的拒绝域为：

 如果满足拒绝域条件，则可以在一定的置信水平上认为线性回归结果显著。

，可以通过查表或者利用软件查出，是检验水平。

当回归效果显著的时候：

可以对系数做置信区间估计，系数的置信区间为

 同理得到相应的置信区间,

点预测和预测区间：

假设是的处的观察结果，是利用回归方程计算的观察值。在给定的置信区间中，

这意味着对于观察值，有的把握认为在区间中。

 例子：合金强度y和其中一个碳含量x有比较密切的关系，今从生产中收集了一批数据如表所示，试拟合一个函数，并且对回归模型进行检验。

x	0.10	0.11	0.12	0.13	0.14	0.15	0.16	0.17	0.18	0.20	0.22	0.24
y	42.0	42.5	45.0	45.5	45.0	47.5	49.0	51.0	50.0	55.0	57.5	59.5

首先我们将图像的散点图画出来，判断使用哪一类方程进行拟合，然后求出方程的未知系数。然后对方程进行拟合优度的检验，最后对系数进行拟合优度的判断。

import pylab as pl
import numpy as np
import statsmodels.api as sm
x = np.array([0.10, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.20, 0.22, 0.24])
y = np.array([42.0, 42.5, 45.0, 45.5, 45.0, 47.5, 49.0, 51.0, 50.0, 55.0, 57.5, 59.5])
def sp(a,b):
    pl.plot(a, b, 'o', color='black')
    pl.show()
def main(d):
    re=sm.formula.ols('y~x',d).fit()
    print(re.summary())
sp(x,y)
x1=x;x1=np.delete(x1,[8,4]);y1=y;y1=np.delete(y1,[8,4])
d={'x':x1,'y':y1}
main(d)

（1）这里首先介绍一下基于公式构建并拟合模型的调用公式，调用格式为：

import statsmodels.api as sm
sm.formula.ols(formula,data=df).fit()

这当中formula是引号括起来的公式，df为数据框或字典格式的数据。使用方法summary获取一个关于此回归的统计检验。

（2）我们在进行数据拟合时通常需要对异常值进行剔除，这里是利用t检验将样本中的异常值进行了剔除：

print(re.outlier_test())

注意：通常使用这个函数，样本数量最好大于20个。

多元线性回归

当自变量的个数多余一个时，回归模型是关于自变量的线性表达形式，则称此模型为多元线性回归，数学模型可以表达为：

 式中：为回归系数；为随机误差，服从正态分布,未知。

回归分析的主要步骤是:

（1）由观测值确定参数（回归系数）的估计值

（2）对线性关系，自变量的显著性进行统计检验

（3）利用回归方程进行预测。

（1） 确定参数

对于模型式中的参数用最小二乘法估计，即选取，使得误差平方和最小。原理和上面一元线性回归模型类似，在这里不过多介绍。

（2） 统计检验

前面是在假设随机变量y和变量j具有线性关系的条件下建立线性回归方程，但是变量y与变量是否为线性关系？所有的变量对变量y是否有影响？都需要做统计检验。

对于总平方,有,是定义的残差平方和，反映随机误差对y的影响；称回归平方和，反映自变量对y的影响。这里是因变量的均值。，其中SST的自由度，SSE的自由度,SSR的自由度。

因变量y与自变量之间是否存在线性规关系是需要检验的，显然如果所有都很小，y和的线性关系就不明显，所以我们可令原假设为

:

当成立时有分解式定义的SSE,SSR满足

~

在显著性水平上，有上分位数,若回归方程效果显著，若,回归效果不显著。

注：

1.当y和的关系不明显时，可能存在非线性关系，例如平方关系。

2.这里的F分布检验直接针对这个模型的所有系数，也就说当这个模型的回归效果显著时，依旧有可能某个自变量与y是不相关的，因此我们还需要对各个回归系数做假设检验，这里使用的是t检验。

上面在一元回归中已经有过定义在这里就不在详述。

（3）回归方程的预测

对于需要预测的自变量，我们直接将其代入我们建立的模型当中，然后就得到了预测值。

我们也可以进行区间估计。记，

,则的置信度为的预测区间可以近似为这是样本量较大时的预测区间。

statsmodels可以使用两种模式求解回归分析模型：一种是基于公式的模型，一种是基于数组的模型。

基于公式的模型：我们在上面的一元线性回归中已经介绍过。

基于数组的模型：

import statsmodels.api as sm
sm.OLS(y,X).fit()

其中y是因变量的观察值向量，而X为自变量观测值矩阵在添加第一列全部为1的增广矩阵。

例题：某企业为了研究车工的平均年龄，平均文化程度和平均产量之间的变化关系，随机抽取了8个班组，测得的数据如下表所示，试求平均产量对平均工龄，平均文化程度的回归方程。并求和时，y的预测值。

	1	2	3	4	5	6	7	8
x1	7.1	6.8	9.2	11.4	8.7	6.6	10.3	10.6
x2	11.1	10.8	12.4	10.9	9.6	9.0	10.5	12.4
x3	15.4	15	22.8	27.8	19.5	13.1	24.9	26.2

import pylab as pl
import numpy as np
import statsmodels.api as sm
a=np.loadtxt('data10_2.txt')
pl.subplot(121)
pl.plot(a[0],a[2],'o',label='x1',color='black')
pl.plot(a[1],a[2],'*',label='x2',color='black')
pl.legend()
d={'x1':a[0], 'x2':a[1],'y':a[2]}
re=sm.formula.ols('y~x1+x2',d).fit()
print(re.summary())
print(re.mse_resid)
yh=re.predict({'x1':[9,10],'x2':[10,9]})
print(yh)
pre=re.get_prediction(d)
df=pre.summary_frame(alpha=0.05)
dfv=df.values;low,upp=dfv[:,4:].T
r=(upp-low)/2
num=np.arange(1,len(a[0])+1)
pl.subplot(122)
pl.errorbar(num,re.resid,r,fmt='o')
pl.show()

这里是基于公式的代码，不做详细介绍。下面附上基于数组的代码，我个人比较推荐基于公式的算法。

import pylab as pl
import numpy as np
import statsmodels.api as sm
a=np.loadtxt('data10_2.txt')
pl.plot(a[0],a[2],'o',label='x1',color='black')
pl.plot(a[1],a[2],'*',label='x2',color='black')
pl.legend()
X=sm.add_constant(a[:2].T)
re=sm.OLS(a[2],X).fit()
print(re.summary())
print(re.mse_resid)
yh=re.predict(np.array([[1,9,10],[1,10,9]]))
print(yh)
pl.show()

（1）这里的函数add_comstant(x)作用是为矩阵x添加一列全为一的作用。

（2）对于基于数组的函数预测和基于公式的有所不同，它的格式和填入OLS函数的格式相同，第一列全为一，第二列为x1，一次类推。

多项式回归

如果我们做出散点图发现y和x呈较明显的二次（或者高次）函数关系，或者用线性关系模拟的效果不太好，就可以使用多项式回归。随机意义下，一元n次多项式回归的数学模型可以表达为

式中：为随机误差，满足,

例子：某厂生产一种电器的销量y和竞争对手的价格和本场的价格有关，我们给出下表，是此商品在10个城市的销售情况，试根据这些数据建立y和和的关系式，并对得到的模型和系数进行检验。并预测若本厂商品在本市售价为160元，那么竞争对手售价为170元。预测这个商品在本市销量。

	120	140	190	130	155	175	125	145	180	150
	100	110	90	150	210	150	250	270	300	250
	102	100	120	77	46	93	26	69	65	85

我们画出散点图，从中可以看出，y和有线性关系，而y和之间的关系难以确定，所以我们多拟合几种模型，去取残差方差最小的这个模型为最佳模型。

线性模型：

纯二次模型：

交叉二次：

纯二次：

import pylab as pl
import numpy as np
import statsmodels.api as sm
a=np.loadtxt('data10_4.txt')
x1=a[0];x2=a[1];y=a[2];
pl.plot(x1,y,'o',label='x1',color='black')
pl.plot(x2,y,'*',label='x2',)
pl.legend()
pl.show()
d={'x1':x1, 'y1':y, 'x2':x2}
resid=[]
idea=['y~x1+x2','y~x1+x2+I(x1**2)','y~x1*x2','y~x1*x2+I(x1**2)+I(x2**2)']
for i in range(0,4):
    re=sm.formula.ols(idea[i],d).fit()
    resid.append(re.mse_resid)
i=resid.index(max(resid))
re=re=sm.formula.ols(idea[i],d).fit()
print('最佳预测值:',re.predict({'x1':170,'x2':160}))

通过选择残差方差这个准则，确定纯二次多项式模型式是拟合的最好的。

逐步回归

实际中影响因变量的因素可能很多，有的可能关联性很强，而有的可能弱一些，人们总希望从中挑选对因变量影响显著的自变量来建立回归模型，下面来讨论多元线性回归模型的情形。

基本思想：为候选自变量的集合，是集合S的一个子集。设中有l个自变量，由和因变量y构造的回归模型的残差平方和为,则模型的残差方差，n为数据样本容量。所选子集应该使尽可能的小。一般来说回归模型所含的自变量越多，那么残差平方和就越小，如果变量中存在对因变量影响很小的量，那么并不会影响，但是会增大导致分母变大，使得残差方差增大。

当自变量较多时，我们有自主选择自变量的方法：

通常我们最常用的是逐步回归法。

前进法：逐步添加变量，每步添加一个，直到没有变量可以引入为止。

第一步：对所有的变量和因变量做一元线性回归，然后得到每个回归模型的F统计量。如果,则将对应变量引入方程中。

第二步：假设我们第一步共引入了m个变量，然后把这m个变量再逐一和剩下的变量建立回归方程，并将其F值和做对比。

第三步：知道所有引入后所有变量的F值均小于。（P为引入变量的个数） 

重复上述步骤，知道添加完所有变量。

后退法：使用与F统计量等价的T检验量

第一步：用所有自变量和因变量做拟合。

第二步：每一步都把未通过t检验的自变量中选择最小一个删除，如果均通过t检验，可以为了使得最小而删除，最小的变量。

第三步：直到所有的自变量均通过t检验。

逐步回归法：

把变量一个一个的引入，每次引入一个变量就对已选入的变量进行逐个检验，当原引入的变量因为后引入的变量而不再显著时，需要进行剔除。引入一个变量或者删除一个变量都要保证模型中的所有变量均为显著的。直到既无显著变量引入，也无不显著的变量从模型里删除。

要求：引入自变量的显著性水平要高于剔除变量的显著性水平。

Logistic回归模型

通常我们会面临二分类问题，比方说我们去医院抽血化验，我们会得到一张化验单上面会标注血液中各成分的含量，医生需要判断我们的身体到底有病或者没病。

这个时候就是一个二分类问题，而Logistic通常给出的是一个概率值，如果有病1，为没病为0的话，我们对一个样本进行预测会得到一个在0~1之间的预测值，如果这个值大于0.5我们就认为他有病，反之则无。

调用格式为：

import statsmodels.api as sm
md=sm.formula.logit('y~x1+x2+x3',d).fit()

调用格式基本和前面的线性回归类似。