ols模型
-
2023年3月1日发(作者：双水解反应方程式)【超详细】多元线性回归模型statsmodels_ols
多元线性模型的主要作⽤：（主要进⾏预测）
通过建模来拟合我们所提供的或是收集到的这些因变量和⾃变量的数据，收集到的数据拟合之后来进⾏参数估计。参数估计的⽬的主要是来
估计出模型的偏回归系数的值。估计出来之后就可以再去收集新的⾃变量的数据去进⾏预测，也称为样本量
importpandasaspd
importnumpyasnp
m#实现了类似于⼆元中的统计模型，⽐如ols普通最⼩⼆乘法
ms#实现了统计⼯具，⽐如t检验、F检验...
mf
importscipy
(991)#随机数种⼦
x1=(0,0.4,100)#⽣成符合正态分布的随机数(均值,标准差,所⽣成随机数的个数)
x2=(0,0.6,100)
x3=(0,0.2,100)
eps=(0,0.05,100)#⽣成噪声数据，保证后⾯模拟所⽣成的因变量的数据⽐较接近实际的环境
X=np.c_[x1,x2,x3]#调⽤c_函数来⽣成⾃变量的数据的矩阵，按照列进⾏⽣成的；100×3的矩阵
beta=[0.1,0.2,0.7]#⽣成模拟数据时候的系数的值
y=(X,beta)+eps#点积+噪声
X_model=_constant(X)#add_constant给矩阵加上⼀列常量1，主要⽬的：便于估计多元线性回归模型的截距，也是便于后⾯进⾏参数估计时的计算
model=(y,X_model)#调⽤OLS普通最⼩⼆乘
results=()#fit拟合，主要功能就是进⾏参数估计，参数估计的主要⽬的是估计出回归系数，根据参数估计结果来计算统计量，这些统计量主要的⽬的
就是对我们模型的有效性或是显著性⽔平来进⾏验证
y()#summary⽅法主要是为了显⽰拟合的结果
⼀、回归系数的参数估计
1.1最⼩⼆乘法实现参数估计——估计⾃变量X的系数
'''
回归系数的计算：X转置乘以X，对点积求逆后，再点乘X转置，最后点乘y
'''
beta_hat=((((X_model.T,X_model)),X_model.T),y)#回归系数
print('回归系数:',(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²
主要作⽤：检验回归模型的显著性
#因变量的回归值=(X_model,系数向量)
y_hat=(X_model,beta_hat)#回归值（拟合值）的计算
y_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(**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²<6说明显著性不⾏了 1.3F检验 F=(ssr/3)/(sse/(100-3-1)); print('F统计量：',round(F,1)) #累积分布函数叫cdf()，调⽤⽅法和下⾯残差函数调⽤⽅法⼀样；注意：cdf+sf算出来的值为1 F_p=(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=#残差 #res=y-y_hat sigma_res=(res)#残差标准差 var_res=(res)#残差⽅差 ll=-(100/2)*(2*)-(100/2)*(var_res)-100/2 print('对数似然：',round(ll,2))#保留两位⼩数 #⾚池信息准则:-2乘以对数似然⽐+2*（参数个数+1）。 #−2ln(L)+2(p+1)（⾚池弘次），其中p为参数个数，ln(L)即ll AIC=-2*ll+2*(3+1) print('AIC：',round(AIC,1)) #贝叶斯信息准则：−2ln(L)+ln(n)∗(p+1),其中ln(L)即llr BIC=-2*ll+(100)*(3+1) print('BIC：',round(BIC,1)) 1.5回归系数标准差 mportt,f C=((X_model.T,X_model))#X倒置点乘X，然后对点集求逆 C_diag=(C)#取出C矩阵对⾓线的值 sigma_unb=(sse/(100-3-1))**(1/2)#残差标准差的⽆偏估计：残差平⽅和/（样本数减参数个数减1） ''' 回归系数标准差stderr的计算： 计算⽅式：残差标准差（⽆偏估计）乘以（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)使⽤残存函数（Survivalfunction）， 等于：1-累积分布函数（1-cdf）； 由于是对t的绝对值进⾏检验，因此需要乘以2。即pt之和 ''' t_const=beta_hat[0]/stdderr_const print('截距项的t值：',round(t_const,3)) p_const=2*(t_const,96) print("P>|t|:",round(p_const,3)) t_x1=beta_hat[1]/stderr_x1 print('x1系数的t值：',round(t_x1,3)) p_t1=2*(t_x1,96) print("P>|t|:",round(p_t1,3)) t_x2=beta_hat[2]/stderr_x2 print('x2系数的t值：',round(t_x2,3)) p_t2=2*(t_x2,96) print("P>|t|:",round(p_t2,3)) t_x3=beta_hat[3]/stderr_x3 print('x3系数的t值：',round(t_x3,3)) p_t3=2*(t_x3,96) print("P>|t|:",round(p_t3,3)) p值>0.05就说明显著性很低，上图截距项的p值=0.329远⼤于0.05，所以不能拒绝原假设，也就是可以去掉截距项 1.7回归系数的置信区间 ''' 回归系数置信区间计算公式： betahat_(1-alpha/2,n-p-1)*sigma_unb*C_diag[0]**(1/2) <=beta_i<= betahat_i+(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来计算分位数的。调⽤进⾏计算， 实际就是cdf函数的逆函数，计算的是百分位数。 ppf和cpf的区别：ppf其实就是通过概率来求随机变量的值；cdf是通过随机变量的值来求概率的 sigma_unb*C_diag[0]**(1/2)就是前⾯计算系数标准差的公式，这⾥可以直接⽤ 系数标准差代替。 特别要注意：要使⽤残差标准差的⽆偏估计，即sigma_unb。 ppf是percentpointfunction的简称，⽤来计算百分位数。 根据置信区间的理论，t分布的概率密度函数关于y轴对称，其均值为0 其显著性⽔平为alpha（置信度为1-alpha）的概率分布函数的分位数，⼀般通过1-alpha/2进⾏计算 ''' const_inter_left=beta_hat[0]-(0.975,96)*sigma_unb*C_diag[0]**(1/2)#置信下限 const_inter_right=beta_hat[0]+(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]-(0.975,96)*sigma_unb*C_diag[1]**(1/2) beta1_inter_right=beta_hat[1]+(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]-(0.975,96)*sigma_unb*C_diag[2]**(1/2) beta2_inter_right=beta_hat[2]+(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]-(0.975,96)*sigma_unb*C_diag[3]**(1/2) beta3_inter_right=beta_hat[3]+(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=((((())**4))/100)/((((())**2)/100)**(1/2)-3#Fisher定义峰度 #Pearson定义峰度 res_kurt=((((())**4))/100)/((((())**2)/100)**2) res_skew=((())**3)/((sigma_res**3)*100) #res_skew= ''' #或者使⽤scipy的函数计算峰度和偏度 mportkurtosis,skew res_kurt=kurtosis(res,fisher=False)#计算峰度，fisher=False表⽰按照pearson定义的⽅式计算峰度 res_skew=skew(res)#计算偏度 #也可以⽤pandas计算偏度和峰度,三种计算⽅式都有⼀定误差 importpandasaspd resFrame=(res) res_kurt=() res_skew=() ''' print('残差峰度(Pearson定义)：',round(res_kurt,3)) print('残差偏度：',round(res_skew,3)) 1.9Jarque-Bera检验与Omnibus检验 1.9.1Jarque-Bera检验 ''' #Jarque-Bera检验，使⽤scipy mportjarque_bera jb_test=jarque_bera(res) ''' '''#使⽤进⾏Jarque-Bera检验 jb_test=_bera(res)#返回值是个命名元组，包含J_B值及其P值 ''' #全⼿动计算J_B值及其P值 mportchi2 #雅克贝拉值=样本量/6*(偏度²+峰度²/4) jb_value=100/6*(res_skew**2+(res_kurt-3)**2/4)#-3是因为⽤的Fisher定义的峰度值，⽽不是⽪尔逊定义的 jb_p=(jb_value,2)#计算雅克贝拉检验其实就是计算雅克贝拉值的分布的残存函数的值，其中的2是卡⽅分布的⾃由度；卡⽅分布的⾃由度是与构成卡⽅ 分布的数据⾥⾯的参数个数有关。上式使⽤的是两个参数：偏度和峰度，因此这⾥⾃由度是2 print('Jarque-Bera(JB):',round(jb_value,3)) print('Prob(JB):',round(jb_p,3)) 1.9.2Omnibus检验 #omnibus检验，使⽤statsmodels omnibus_test=_normtest(res)#omnibus检验 print('Omnibus：',round(omnibus_tic,3))#omnibus的统计量 print('Prob(Omnibus):',round(omnibus_,3))#omnibus的p值 ''' Omnibus检验的具体步骤： （1）计算残差的偏度检验值和峰度检验值。 （2）求出⼆者平⽅和。 （3）以平⽅和和⾃由度2为参数调⽤卡⽅分布的残存函数。 （4）平⽅和是Omnibus统计量的值，残存函数返回值是Omnibus统计量的P值。 从此⽰例可以看出scipy科学与统计计算功能之强⼤！ ''' mportnormaltest,skewtest,kurtosistest,skew,kurtosis,chi2 #normaltest(res)#此函数直接进⾏Omnibus检验 s,_=skewtest(res)#注意：偏度检验和偏度并不是⼀回事 k,_=kurtosistest(res)#峰度检验和峰度也不是⼀回事 k2=s*s+k*k print('------------------⼿动计算Omnibus检验----------------------') print('Omnibus:',(k2,3)) #print('Prob(Omnibus):',((k2,2),3))#通过卡⽅分布的残存函数计算Omnibus的P值。 1.10Dubin-Watson检验 ⾏与⾏之间⾃相关，2附近表⽰没有相关 若越接近于0或越接近于4都表⽰存在相关 ''' Durbin-Watson检验：越接近2，表⽰残差越接近正态分布 #直接调⽤_watson()函数 dw=_watson(res) ''' #原始计算公式：残差的差值平⽅和除以残差平⽅和 diff_resids=(res) dw=(diff_resids**2)/(res**2) print('Durbin-Watson:',round(dw,3)) 1.11条件数. 900以内为佳，>900说明有共线性 ''' 条件数（.）的计算步骤如下： （1）获取增加常量1向量后的⾃变量矩阵，X_model （2）计算X_model转置与其本⾝的点积C （3）计算点积的特征值 （4）最⼤特征值/最⼩特征值，然后将结果开平⽅ 条件数是衡量矩阵病态的⼀个指标，理论上该值越⼩越好。 ''' #C=(X_model.T,X_model) eigs=(C)[0]#linalg线性代数的缩写，eigh特征值单词的缩写。eig也可以计算特征值，主要可以计算⼀般的矩阵；eigh计算的矩阵是对称矩阵， 计算出来的特征值会从⼩到⼤进⾏排序 cond=(eigs[-1]/eigs[0])#(最⼤特征值/最⼩特征值)² print('通过特征值计算.:',round(cond,3)) ''' 条件数的另⼀种计算法：sqrt(||C||*||inv(C)||)， 即⽤矩阵C的范数乘以C逆的范数，然后再开平⽅。 ''' cond=((C,ord=2)*((C),ord=2))#注意设置ord=2 print('通过矩阵计算.:',round(cond,3)) 1.12补充：VIF（只适⽤于⼩数据） ''' 多元回归分析的⽅差膨胀系数（VIF）的计算，这⾥直接调⽤StatsModels的函数。 某些教科书认为VIF>10即判定⾃变量存在多重共线性；StatsModels将阈值设为5。 ''' rs_influenceimportvariance_inflation_factor foriinrange([1]): print('x%d的VIF：%0.4f'%(i,variance_inflation_factor(X,i))) #计算第i个⾃变量的VIF，⽐如i=1 k_vars=[1]#提取列数，即⾃变量个数 x_1=X[:,1] mask=(k_vars)!=1 x_not1=X[:,mask] #以第i个⾃变量作为因变量，其他⾃变量作为⾃变量调⽤OLS，然后提取Rsquared_i r_squared_1=(x_1,x_not1).fit().rsquared vif_i1=1./(1.-r_squared_1) print('⼿动计算⾃变量VIF（x1）：',vif_i1) 1.13附录 ''' 说明：查阅有关通过范数计算条件数的资料，都是⽤矩阵的范数乘以该矩阵逆的范数，都⽆开平⽅这个计算步骤。 statsmodels的开平⽅这个计算步骤不知是何⽤意？暂时还没弄明⽩。 ⽰例：/qq_18343569/article/details/50404989 ''' m=([[3,5,0],[2,10,4],[3,4,5]]) #m=ix(a) (m,ord=2)*((m),ord=2)#norm(需要计算的矩阵,ord=2)：计算矩阵的范数,ord表⾯使⽤的谱范数中的2-范数 ''' 泛函分析中的范数计算. 谱范数的2-范数（注意与l2范数区别）：等于矩阵的共轭转置和⾃⾝的点积⽣成矩阵的特征值中的最⼤值开平⽅。 很显然：(matrix,ord=2)等中的ord=2表⽰计算矩阵的2-范数。 ''' eigs_norm=(((m.T,m))[0]) print('矩阵m的2-范数（通过公式计算）：',eigs_()) norm_2=(m,ord=2) print('矩阵m的2-范数（调⽤numpy函数）：',norm_2) ''' l2范数计算与谱范数的2-范数计算有区别：前者表⽰矩阵元素平⽅和，然后再开平⽅，如下例所⽰： ''' l2norm_byhand=((m**2)) print('⼿⼯计算l2norm：',l2norm_byhand) l2norm_bynormfunc=(m)#如果ord为空，缺省即l2范数 print('通过norm函数计算l2norm：',l2norm_bynormfunc) #多维数组，*乘法是逐个元素相乘 m*m m@m#矩阵点积 m1=ix(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²系数、残差