多元线性模型的主要作用：（主要进行预测）

通过建模来拟合我们所提供的或是收集到的这些因变量和自变量的数据，收集到的数据拟合之后来进行参数估计。参数估计的目的主要是来估计出模型的偏回归系数的值。估计出来之后就可以再去收集新的自变量的数据去进行预测，也称为样本量

import pandas as pdimport numpy as npimport statsmodels.api as sm#实现了类似于二元中的统计模型，比如ols普通最小二乘法import statsmodels.stats.api as sms#实现了统计工具，比如t检验、F检验...import statsmodels.formula.api as smfimport scipynp.random.seed(991)#随机数种子x1 = np.random.normal(0,0.4,100)#生成符合正态分布的随机数(均值,标准差,所生成随机数的个数)x2 = np.random.normal(0,0.6,100)x3 = np.random.normal(0,0.2,100)eps = np.random.normal(0,0.05,100)#生成噪声数据，保证后面模拟所生成的因变量的数据比较接近实际的环境X = np.c_[x1,x2,x3]#调用c_函数来生成自变量的数据的矩阵，按照列进行生成的；100×3的矩阵beta = [0.1,0.2,0.7]#生成模拟数据时候的系数的值y = np.dot(X,beta) + eps#点积+噪声X_model = sm.add_constant(X)#add_constant给矩阵加上一列常量1，主要目的：便于估计多元线性回归模型的截距，也是便于后面进行参数估计时的计算model = sm.OLS(y,X_model)#调用OLS普通最小二乘results = model.fit()#fit拟合，主要功能就是进行参数估计，参数估计的主要目的是估计出回归系数，根据参数估计结果来计算统计量，这些统计量主要的目的就是对我们模型的有效性或是显著性水平来进行验证results.summary()#summary方法主要是为了显示拟合的结果

一、回归系数的参数估计

1.1 最小二乘法实现参数估计——估计自变量X的系数

'''回归系数的计算：X转置乘以X，对点积求逆后，再点乘X转置，最后点乘y '''beta_hat = np.dot(np.dot(np.linalg.inv(np.dot(X_model.T,X_model)),X_model.T),y)#回归系数print('回归系数:',np.round(beta_hat,4))#四舍五入取小数点后4位print('回归方程：Y_hat=%0.4f+%0.4f*X1+%0.4f*X2+%0.4f*X3' % (beta_hat[0],beta_hat[1],beta_hat[2],beta_hat[3]))

1.2 决定系数：R²与调整后R²

主要作用：检验回归模型的显著性

#因变量的回归值=np.dot(X_model,系数向量)y_hat = np.dot(X_model,beta_hat)#回归值（拟合值）的计算y_mean = np.mean(y)#求因变量的平均值sst = sum((y-y_mean)**2)#总平方和：即y减去y均值后差的平方和ssr = sum((y_hat-y_mean)**2)#回归平方和:y回归值减去y均值后差的平方和#sse = sum((y-y_hat)**2)#残差平方和:y值减去y回归值之差的平方和sse = sum(results.resid**2)#结果和上面注释了的式子一样，或许有小数点的误差，但基本上可忽略不计R_squared =1 - sse/sst#R²：1减去残差平方和除以总平方和商的差；求解方法二：R²=ssr/sst#按照线性回归模型原理来说：[残差平方和+回归平方和=总平方和]→[R²=ssr/sst]print('R方:',round(R_squared,3))#调整后平方和：100表示样本数据总数(n)，3表示自变量个数(p)adjR_squared =1- (sse/(100-3-1))/(sst/(100-1))#1-(残差的平方和/残差的自由度)/(总平方和/无偏估计)#残差的自由度也是残差方差的无偏估计print('调整后R方:',round(adjR_squared,3))

0≤R²≤1

R²>8说明模型显著性强

6<R²≤8说明模型还行

R²<6说明显著性不行了

1.3 F检验

F = (ssr/3)/(sse/(100-3-1));print('F统计量：',round(F,1))#累积分布函数叫cdf()，调用方法和下面残差函数调用方法一样；注意：cdf+sf算出来的值为1F_p = scipy.stats.f.sf(F,3,96)#使用F分布的残存函数计算P值；sf是残存函数英语单词的缩写；3和96分别是两个自由度print('F统计量的P值：', F_p)

1.4 对数似然、AIC与BIC

#对数似然值计算公式： L=-(n/2)*ln(2*pi)-(n/2)*ln(sse/n)-n/2；sse/n就是方差res = results.resid#残差#res = y-y_hatsigma_res = np.std(res) #残差标准差var_res = np.var(res) #残差方差ll = -(100/2)*np.log(2*np.pi)-(100/2)*np.log(var_res)-100/2print('对数似然：', round(ll,2))#保留两位小数#赤池信息准则:-2乘以对数似然比+2*（参数个数+1）。#−2ln(L)+2(p+1)（赤池弘次），其中p为参数个数，ln(L)即llAIC = -2*ll + 2*(3+1)print('AIC：',round(AIC,1))# 贝叶斯信息准则：−2ln(L)+ln(n)∗(p+1),其中ln(L)即llrBIC = -2*ll+np.log(100)*(3+1) print('BIC：',round(BIC,1))

1.5 回归系数标准差

from scipy.stats import t,fC = np.linalg.inv(np.dot(X_model.T,X_model))#X倒置点乘X，然后对点集求逆C_diag = np.diag(C)#取出C矩阵对角线的值sigma_unb= (sse/(100-3-1))**(1/2)#残差标准差的无偏估计：残差平方和/（样本数减参数个数减1）'''回归系数标准差std err的计算：计算方式：残差标准差（无偏估计）乘以（C矩阵对角线上对应值的平方根）'''stdderr_const = sigma_unb*(C_diag[0]**(1/2))#常数项（截距）的标准差，对应C_diag[0]print('常数项（截距）的标准差：',round(stdderr_const,3))stderr_x1 = sigma_unb*(C_diag[1]**(1/2))#第一个系数对应C_diag[1]print('beta1的标准差：',round(stderr_x1,3))stderr_x2 = sigma_unb*(C_diag[2]**(1/2))#第二个系数对应C_diag[2]print('beta2的标准差：',round(stderr_x2,3))stderr_x3 = sigma_unb*(C_diag[3]**(1/2))#第三个系数对应C_diag[3]print('beta3的标准差：',round(stderr_x3,3))

print('C矩阵：\n', C)print('\nC矩阵的对角线元素：',C_diag)

1.6 回归系数的显著性检验（t检验）

t值足够小就可以认为回归方程的系数是具有显著性的，显著性水平是比较高的，否则就可以认为这个回归系数的显著性不高

'''回归方程的显著性检验(1) t检验：beta_hat[i](相应的回归系数)除以相应系数标准差(2) 使用scipy.stats.t.sf残存函数（Survival function），等于：1-累积分布函数（1-cdf）；由于是对t的绝对值进行检验，因此需要乘以2。即p<-t与p>t之和 '''t_const = beta_hat[0]/stdderr_constprint('截距项的t值：',round(t_const,3))p_const = 2*t.sf(t_const,96)print("P>|t|:",round(p_const,3))t_x1 = beta_hat[1]/stderr_x1print('x1系数的t值：',round(t_x1,3))p_t1 = 2*t.sf(t_x1,96)print("P>|t|:",round(p_t1,3))t_x2 = beta_hat[2]/stderr_x2print('x2系数的t值：',round(t_x2,3))p_t2 = 2*t.sf(t_x2,96)print("P>|t|:",round(p_t2,3))t_x3 = beta_hat[3]/stderr_x3print('x3系数的t值：',round(t_x3,3))p_t3 = 2*t.sf(t_x3,96)print("P>|t|:",round(p_t3,3))

p值>0.05就说明显著性很低，上图截距项的p值=0.329远大于0.05，所以不能拒绝原假设，也就是可以去掉截距项

1.7 回归系数的置信区间

'''回归系数置信区间计算公式：betahat_i - t.ppf(1-alpha/2,n-p-1)*sigma_unb*C_diag[0]**(1/2)< =beta_i <=betahat_i + t.ppf(1-alpha/2,n-p-1)*sigma_unb*C_diag[0]**(1/2)t(1-alpha/2,n-p-1)是t值在0.025，自由度为n-p的分位数，由于下分位数是负值，所以这里使用0.975来计算分位数的。调用scipy.stats.t.ppf进行计算，实际就是cdf函数的逆函数，计算的是百分位数。ppf和cpf的区别：ppf其实就是通过概率来求随机变量的值；cdf是通过随机变量的值来求概率的sigma_unb*C_diag[0]**(1/2)就是前面计算系数标准差的公式，这里可以直接用系数标准差代替。特别要注意：要使用残差标准差的无偏估计，即sigma_unb。ppf是percent point function的简称，用来计算百分位数。根据置信区间的理论，t分布的概率密度函数关于y轴对称，其均值为0其显著性水平为alpha（置信度为1-alpha）的概率分布函数的分位数，一般通过1-alpha/2进行计算'''const_inter_left = beta_hat[0] - t.ppf(0.975,96)*sigma_unb*C_diag[0]**(1/2)#置信下限const_inter_right = beta_hat[0] + t.ppf(0.975,96)*sigma_unb*C_diag[0]**(1/2)#置信上限print('常数项的置信区间是: [%0.3f,%0.3f]'%(const_inter_left,const_inter_right))beta1_inter_left = beta_hat[1] - t.ppf(0.975,96)*sigma_unb*C_diag[1]**(1/2)beta1_inter_right = beta_hat[1] + t.ppf(0.975,96)*sigma_unb*C_diag[1]**(1/2)print('x1系数的置信区间是: [%0.3f,%0.3f]'%(beta1_inter_left,beta1_inter_right))beta2_inter_left = beta_hat[2] - t.ppf(0.975,96)*sigma_unb*C_diag[2]**(1/2)beta2_inter_right = beta_hat[2] + t.ppf(0.975,96)*sigma_unb*C_diag[2]**(1/2)print('x2系数的置信区间是: [%0.3f,%0.3f]'%(beta2_inter_left,beta2_inter_right))beta3_inter_left = beta_hat[3] - t.ppf(0.975,96)*sigma_unb*C_diag[3]**(1/2)beta3_inter_right = beta_hat[3] + t.ppf(0.975,96)*sigma_unb*C_diag[3]**(1/2)print('x3系数的置信区间是: [%0.3f,%0.3f]' % (beta3_inter_left,beta3_inter_right))

1.8 峰度与偏度

'''1、注意峰度的定义方式有两种：一是Fisher定义，正态分布值为0；另一个是Pearson定义，正态分布值为3。StatsModels使用的是Pearson定义。按照Fisher定义，逢峰度=0表示正好符合正正态分布，大于0表示峰比较尖，反之表示比较平。2、对于偏度，偏度值大于0则为正偏态或左偏态；小于零则表示负偏态或右偏态。'''#kurtosis = ((np.sum((res-res.mean())**4))/100)/((np.sum ((res-res.mean())**2)/100)**(1/2)-3 #Fisher定义峰度#Pearson定义峰度res_kurt = ((np.sum((res-res.mean())**4))/100)/((np.sum ((res-res.mean())**2)/100)**2)res_skew = np.sum((res-res.mean())**3)/((sigma_res**3)*100)#res_skew = '''#或者使用scipy的函数计算峰度和偏度from scipy.stats import kurtosis,skewres_kurt = kurtosis(res,fisher=False) #计算峰度，fisher=False表示按照pearson定义的方式计算峰度res_skew = skew(res) #计算偏度#也可以用pandas计算偏度和峰度,三种计算方式都有一定误差import pandas as pdresFrame = pd.Series(res)res_kurt = resFrame.kurt()res_skew = resFrame.skew()'''print('残差峰度(Pearson定义)：',round(res_kurt,3))print('残差偏度：',round(res_skew,3))

1.9 Jarque-Bera检验与Omnibus检验

1.9.1 Jarque-Bera检验

'''#Jarque-Bera检验，使用scipyfrom scipy.stats import jarque_berajb_test = jarque_bera(res)''''''#使用statsmodels.stats.api进行Jarque-Bera检验jb_test = sms.jarque_bera(res)#返回值是个命名元组，包含J_B值及其P值'''#全手动计算J_B值及其P值from scipy.stats import chi2#雅克贝拉值=样本量/6*(偏度²+峰度²/4)jb_value = 100 / 6 * (res_skew**2 + (res_kurt - 3)**2 / 4)#-3是因为用的Fisher定义的峰度值，而不是皮尔逊定义的jb_p = chi2.sf(jb_value, 2)#计算雅克贝拉检验其实就是计算雅克贝拉值的分布的残存函数的值，其中的2是卡方分布的自由度；卡方分布的自由度是与构成卡方分布的数据里面的参数个数有关。上式使用的是两个参数：偏度和峰度，因此这里自由度是2print('Jarque-Bera (JB): ',round(jb_value,3))print('Prob(JB): ',round(jb_p,3))

1.9.2 Omnibus检验

#omnibus检验，使用statsmodelsomnibus_test = sms.omni_normtest(res) #omnibus检验print('Omnibus：',round(omnibus_test.statistic,3))#omnibus的统计量print('Prob(Omnibus): ',round(omnibus_test.pvalue,3))#omnibus的p值'''Omnibus检验的具体步骤：（1）计算残差的偏度检验值和峰度检验值。（2）求出二者平方和。（3）以平方和和自由度2为参数调用卡方分布的残存函数。（4）平方和是Omnibus统计量的值，残存函数返回值是Omnibus统计量的P值。从此示例可以看出scipy科学与统计计算功能之强大！'''from scipy.stats import normaltest,skewtest,kurtosistest,skew,kurtosis,chi2#normaltest(res)#此函数直接进行Omnibus检验s, _ = skewtest(res)#注意：偏度检验和偏度并不是一回事k, _ = kurtosistest(res)#峰度检验和峰度也不是一回事k2 = s*s + k*kprint('------------------手动计算Omnibus检验----------------------')print('Omnibus: ', np.round(k2,3))#print('Prob(Omnibus):', np.round(chi2.sf(k2,2),3))#通过卡方分布的残存函数计算Omnibus的P值。

1.10 Dubin-Watson检验

行与行之间自相关，2附近表示没有相关

若越接近于0或越接近于4都表示存在相关

'''Durbin-Watson检验：越接近2，表示残差越接近正态分布#直接调用sms.durbin_watson()函数dw = sms.durbin_watson(res)'''#原始计算公式：残差的差值平方和除以残差平方和diff_resids = np.diff(res)dw = np.sum(diff_resids**2) / np.sum(res**2)print('Durbin-Watson: ',round(dw,3))

1.11 条件数Cond. No.

900以内为佳，>900说明有共线性

'''条件数（Cond. No.）的计算步骤如下：（1）获取增加常量1向量后的自变量矩阵，X_model（2）计算X_model转置与其本身的点积C（3）计算点积的特征值（4）最大特征值/最小特征值，然后将结果开平方条件数是衡量矩阵病态的一个指标，理论上该值越小越好。'''#C= np.dot(X_model.T,X_model)eigs = np.linalg.eigh(C)[0] #linalg线性代数的缩写，eigh特征值单词的缩写。eig也可以计算特征值，主要可以计算一般的矩阵；eigh计算的矩阵是对称矩阵，计算出来的特征值会从小到大进行排序cond = np.sqrt(eigs[-1]/eigs[0])#(最大特征值/最小特征值)²print('通过特征值计算Cond. No.: ',round(cond,3))'''条件数的另一种计算法：sqrt(||C||*||inv(C)||)，即用矩阵C的范数乘以C逆的范数，然后再开平方。'''cond = np.sqrt(np.linalg.norm(C,ord=2)*np.linalg.norm(np.linalg.inv(C),ord=2))#注意设置ord=2print('通过矩阵计算Cond. No.: ',round(cond,3))

1.12 补充：VIF（只适用于小数据）

'''多元回归分析的方差膨胀系数（VIF）的计算，这里直接调用StatsModels的函数。某些教科书认为VIF>10即判定自变量存在多重共线性；StatsModels将阈值设为5。'''from statsmodels.stats.outliers_influence import variance_inflation_factorfor i in range(X.shape[1]):print('x%d的VIF：%0.4f'%(i,variance_inflation_factor(X,i)))#计算第i个自变量的VIF，比如i=1k_vars = X.shape[1]#提取列数，即自变量个数x_1 = X[:, 1]mask = np.arange(k_vars) != 1x_not1 = X[:, mask]#以第i个自变量作为因变量，其他自变量作为自变量调用OLS，然后提取Rsquared_ir_squared_1 = sm.OLS(x_1, x_not1).fit().rsquaredvif_i1= 1. / (1. - r_squared_1)print('手动计算自变量VIF（x1）：',vif_i1)

1.13 附录

'''说明：查阅有关通过范数计算条件数的资料，都是用矩阵的范数乘以该矩阵逆的范数，都无开平方这个计算步骤。statsmodels的开平方这个计算步骤不知是何用意？暂时还没弄明白。示例：/qq_18343569/article/details/50404989'''m=np.array([[3,5,0],[2,10,4],[3,4,5]])#m = np.asmatrix(a)np.linalg.norm(m,ord=2)*np.linalg.norm(np.linalg.inv(m),ord=2)#norm(需要计算的矩阵,ord=2)：计算矩阵的范数,ord表面使用的谱范数中的2-范数

'''泛函分析中的范数计算.谱范数的2-范数（注意与l2范数区别）：等于矩阵的共轭转置和自身的点积生成矩阵的特征值中的最大值开平方。很显然：np.linalg.norm(matrix,ord=2)等中的ord=2表示计算矩阵的2-范数。'''eigs_norm = np.sqrt(np.linalg.eigh(np.dot(m.T,m))[0])print('矩阵m的2-范数（通过公式计算）：',eigs_norm.max())norm_2 = np.linalg.norm(m,ord=2)print('矩阵m的2-范数（调用numpy函数）：',norm_2)

'''l2范数计算与谱范数的2-范数计算有区别：前者表示矩阵元素平方和，然后再开平方，如下例所示：'''l2norm_byhand = np.sqrt(np.sum(m**2))print('手工计算l2norm：',l2norm_byhand)l2norm_bynormfunc = np.linalg.norm(m)#如果ord为空，缺省即l2范数print('通过norm函数计算l2norm：',l2norm_bynormfunc)

#多维数组，*乘法是逐个元素相乘m*m

m@m#矩阵点积

m1 = np.asmatrix(m)#如果将多维数组转换成矩阵，*乘法就是矩阵乘法——点积m1*m1

#@也是矩阵乘法，对于多维数组和矩阵是一样的。m1@m1

#多维数组支持转置但是不支持共轭转置。对于实数矩阵，转置和共轭转置的结果相同print('多维数组的转置：', m.T)m.H #求数组的共轭转置的话就会发生异常

#矩阵支持共轭转置m1.H

多元线性回归分析：归因和预测问题，y连续，x分类glm，R²

1、图形：散点图（y和x之间的关系）异常、线性关系、稀疏问题、相关强弱

2、相关：系数大小：0.1以上才有关系、

3、回归：R²(0.35-0.5)，系数p，共线性超过900就要处理

4、残差：正态、异方差、异常值、内生性。

5、判断什么是主要原因1？——内生

6、处理内生性，——R²系数，残差

7、判断什么是主要原因2？——异常

8、处理异常值，——R²系数、残差