python3 线性回归检验2

#====图示法完成方差齐性的判断=====
#标准化残差与预测值之间的散点图
plt.scatter(fit2.predict(),(fit2.resid-fit2.resid.mean())/fit2.resid.std())
plt.xlabel('预测值')
plt.ylabel('标准化残差')
#添加水平参考线
plt.axhline(y = 0,color = 'r',linewidth = 2)
plt.show()
#====统计法完成方差齐性的判断=====
#White`s Test
sm.stats.diagnostic.het_white(fit2.resid,exog = fit2.model.exog)
#Breush-Pagan
sm.stats.diagnostic.het_breushpagan(fit2.resid,exog_het = fit2.model.exog)

#==========残差与x的关系=======
plt.subplot(231)
plt.scatter(ccpp_outliers.AT,(fit2.resid-fit2.resid.mean())/fit2.resid.std())
plt.xlabel('AT')
plt.ylabel('标准化残差')
plt.axhline(color = 'red',linewidth = 2)

plt.subplot(232)
plt.scatter(ccpp_outliers.V,(fit2.resid-fit2.resid.mean())/fit2.resid.std())
plt.xlabel('V')
plt.ylabel('标准化残差')
plt.axhline(color = 'red',linewidth = 2)

plt.subplot(233)
plt.scatter(ccpp_outliers.AP,(fit2.resid-fit2.resid.mean())/fit2.resid.std())
plt.xlabel('AP')
plt.ylabel('标准化残差')
plt.axhline(color = 'red',linewidth = 2)

plt.subplot(234)
plt.scatter(np.power(ccpp_outliers.AT,2),(fit2.resid-fit2.resid.mean())/fit2.resid.std())
plt.xlabel('AT^2')
plt.ylabel('标准化残差')
plt.axhline(color = 'red',linewidth = 2)

plt.subplot(235)
plt.scatter(np.power(ccpp_outliers.V,2),(fit2.resid-fit2.resid.mean())/fit2.resid.std())
plt.xlabel('V^2')
plt.ylabel('标准化残差')
plt.axhline(color = 'red',linewidth = 2)

plt.subplot(236)
plt.scatter(np.power(ccpp_outliers.AP,2),(fit2.resid-fit2.resid.mean())/fit2.resid.std())
plt.xlabel('AP^2')
plt.ylabel('标准化残差')
plt.axhline(color = 'red',linewidth = 2)

#设置子图之间的水平间距和高度间距
plt.subplots_adjust(hspace = 0.3,wspace = 0.3)
plt.show()

# ====== 统计法完成方差齐性的判断 ======
# White's Test
sm.stats.diagnostic.het_white(fit2.resid, exog = fit2.model.exog)
# Breusch-Pagan
sm.stats.diagnostic.het_breushpagan(fit2.resid, exog_het = fit2.model.exog)
# 三种权重
w1 = 1/np.abs(fit2.resid)
w2 = 1/fit2.resid**2
ccpp_outliers['loge2'] = np.log(fit2.resid**2)
# 第三种权重
model = sm.formula.ols('loge2~AT+V+AP', data = ccpp_outliers).fit()
w3 = 1/(np.exp(model.predict()))

# 建模,此处重点注意weights要写对！！！！
fit3 = sm.formula.wls('PE~AT+V+AP', data = ccpp_outliers, weights = w1).fit()
# 异方差检验
het3 = sm.stats.diagnostic.het_breushpagan(fit3.resid, exog_het = fit3.model.exog)
# AIC
fit3.aic

fit4 = sm.formula.wls('PE~AT+V+AP', data = ccpp_outliers, weights = w2).fit()
het4 = sm.stats.diagnostic.het_breushpagan(fit4.resid, exog_het = fit4.model.exog)
fit4.aic

fit5 = sm.formula.wls('PE~AT+V+AP', data = ccpp_outliers, weights = w3).fit()
het5 = sm.stats.diagnostic.het_breushpagan(fit5.resid, exog_het = fit5.model.exog)
fit5.aic

# fit2模型
het2 = sm.stats.diagnostic.het_breushpagan(fit2.resid, exog_het = fit2.model.exog)
fit2.aic

print('fit2模型异方差检验统计量：%.2f，P值为%.4f：' %(het2[0],het2[1]))
print('fit3模型异方差检验统计量：%.2f，P值为%.4f：' %(het3[0],het3[1]))
print('fit4模型异方差检验统计量：%.2f，P值为%.4f：' %(het4[0],het4[1]))
print('fit5模型异方差检验统计量：%.2f，P值为%.4f：\n' %(het5[0],het5[1]))

print('fit2模型的AIC：%.2f' %fit2.aic)
print('fit3模型的AIC：%.2f' %fit3.aic)
print('fit4模型的AIC：%.2f' %fit4.aic)
print('fit5模型的AIC：%.2f' %fit5.aic)