多元线性模型的主要作用：（主要进行预测）
通过建模来拟合我们所提供的或是收集到的这些因变量和自变量的数据，收集到的数据拟合之后来进行参数估计。参数估计的目的主要是来估计出模型的偏回归系数的值。估计出来之后就可以再去收集新的自变量的数据去进行预测，也称为样本量
import pandas as pd
import numpy as np
import statsmodels.api as sm#实现了类似于二元中的统计模型，比如ols普通最小二乘法
import statsmodels.stats.api as sms#实现了统计工具，比如t检验、F检验...
import statsmodels.formula.api as smf
import scipy
np.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算出来的值为1
F_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_hat
sigma_res = np.std(res) #残差标准差
var_res = np.var(res) #残差方差
ll = -(100/2)*np.log(2*np.pi)-(100/2)*np.log(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+np.log(100)*(3+1)
print('BIC：',round(BIC,1))
1.5 回归系数标准差
from scipy.stats import t,f
C = 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_const
print('截距项的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_x1
print('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_x2
print('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_x3
print('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,skew
res_kurt = kurtosis(res,fisher=False) #计算峰度，fisher=False表示按照pearson定义的方式计算峰度
res_skew = skew(res) #计算偏度
#也可以用pandas计算偏度和峰度,三种计算方式都有一定误差
import pandas as pd
resFrame = 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检验，使用scipy
from scipy.stats import jarque_bera
jb_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是卡方分布的自由度；卡方分布的自由度是与构成卡方分布的数据里面的参数个数有关。上式使用的是两个参数：偏度和峰度，因此这里自由度是2
print('Jarque-Bera (JB): ',round(jb_value,3))
print('Prob(JB): ',round(jb_p,3))
1.9.2 Omnibus检验
#omnibus检验，使用statsmodels
omnibus_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*k
print('------------------手动计算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=2
print('通过矩阵计算Cond. No.: ',round(cond,3))
1.12 补充：VIF（只适用于小数据）
'''
多元回归分析的方差膨胀系数（VIF）的计算，这里直接调用StatsModels的函数。
某些教科书认为VIF>10即判定自变量存在多重共线性；StatsModels将阈值设为5。
'''
from statsmodels.stats.outliers_influence import variance_inflation_factor
for i in range(X.shape[1]):
print('x%d的VIF：%0.4f'%(i,variance_inflation_factor(X,i)))
#计算第i个自变量的VIF，比如i=1
k_vars = X.shape[1]#提取列数，即自变量个数
x_1 = X[:, 1]
mask = np.arange(k_vars) != 1
x_not1 = X[:, mask]
#以第i个自变量作为因变量，其他自变量作为自变量调用OLS，然后提取Rsquared_i
r_squared_1 = sm.OLS(x_1, x_not1).fit().rsquared
vif_i1= 1. / (1. - r_squared_1)
print('手动计算自变量VIF（x1）：',vif_i1)
1.13 附录
'''
说明：查阅有关通过范数计算条件数的资料，都是用矩阵的范数乘以该矩阵逆的范数，都无开平方这个计算步骤。
statsmodels的开平方这个计算步骤不知是何用意？暂时还没弄明白。
示例：https://blog.csdn.net/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²系数、残差
扫一扫
814
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
