python cox回归_python3 线性回归验证方法

#-*- coding: utf-8 -*-

import pandas as pd

import numpy as np

from patsy.highlevel import dmatrices

#2.7里面是from patsy import dmatrices

from statsmodels.stats.outliers_influence import variance_inflation_factor

import statsmodels.api as sm

import scipy.stats as stats

from sklearn.metrics import mean_squared_error

import seaborn as sns

import matplotlib.pyplot as plt

import matplotlib.mlab as mlab

import matplotlib

#数据获取

ccpp = pd.read_excel('CCPP.xlsx')

ccpp.describe()

#绘制各变量之间的散点图

sns.pairplot(ccpp)

plt.show()

#发电量(PE)与自变量之间的相关系数

a = ccpp.corrwith(ccpp.PE)

print(a)

#将因变量PE，自变量AT,V,AP和截距项(值为1的1维数值)以数据框的形式组合起来

y,x = dmatrices('PE~AT+V+AP',data = ccpp,return_type = 'dataframe')

#构造空的数据框

vif = pd.DataFrame()

vif[""VIF Factor""] = [variance_inflation_factor(x.values,i) for i in range(x.shape[1])]

vif[""features""] = x.columns

print (vif)

#构建PE与AT,V和AP之间的线性模型

fit = sm.formula.ols('PE~AT+V+AP',data=ccpp).fit()

b = fit.summary()

# print(b)

#计算模型的RMSE值

pred = fit.predict()

c = np.sqrt(mean_squared_error(ccpp.PE,pred))

print(c)

#离群点检验

outliers = fit.get_influence()

#高杠杆值点(帽子矩阵)

leverage = outliers.hat_matrix_diag

#dffits值

dffits = outliers.dffits[0]

#学生化残差

resid_stu = outliers.resid_studentized_external

#cook距离

cook = outliers.cooks_distance[0]

#covratio值

covratio = outliers.cov_ratio

#将上面的几种异常值检验统计量与原始数据集合并

contat1 = pd.concat([pd.Series(leverage,name = 'leverage'),pd.Series(dffits,name ='dffits'),

pd.Series(resid_stu,name = 'resid_stu'),pd.Series(cook,name = 'cook'),

pd.Series(covratio,name ='covratio'),],axis = 1)

ccpp_outliers = pd.concat([ccpp,contat1],axis = 1)

d = ccpp_outliers.head()

print(d)

#计算异常值数量的比例

outliers_ratio = sum(np.where((np.abs(ccpp_outliers.resid_stu)>2),1,0))/ccpp_outliers.shape[0]

e = outliers_ratio

print(e)

#删除异常值

ccpp_outliers = ccpp_outliers.loc[np.abs(ccpp_outliers.resid_stu)<=2,]

#重新建模

fit2 = sm.formula.ols('PE~AT+V+AP',data = ccpp_outliers).fit()

f = fit2.summary()

# print(f)

pred2 = fit2.predict()

g = np.sqrt(mean_squared_error(ccpp_outliers.PE,pred2))

print(g)

#

#残差的正态性检验(直方图法)

resid = fit2.resid

#中文和负号的正常显示

# plt.rcParams['font.sans=serif'] = ['Microsoft YaHei']

plt.rcParams['font.sans-serif'] = ['SimHei']

# plt.rcParams['font.sans=serif'] = 'sans-serif'

plt.rcParams['axes.unicode_minus'] = False

plt.hist(resid,bins = 100,normed = True,color = 'steelblue',edgecolor = 'k')

#设置坐标轴标签和标题

plt.title('残差直方图')

plt.ylabel('密度值')

#生成正态曲线的数据

x1 = np.linspace(resid.min(),resid.max(),1000)

normal = mlab.normpdf(x1,resid.mean(),resid.std())

#绘制正态分布曲线

plt.plot(x1,normal,'r-',linewidth = 2,label = '正态分布曲线')

#生成核密度曲线的数据

kde = mlab.GaussianKDE(resid)

x2 = np.linspace(resid.min(),resid.max(),1000)

#绘制核密度曲线

plt.plot(x2,kde(x2),'k-',linewidth = 2,label = '核密度曲线')

#去除图形顶部边界和右边界的刻度

plt.tick_params(top = 'off',right = 'off')

#显示图例

plt.legend(loc='best')

#显示图形

plt.show()

#生成的正态曲线的数据

pp_qq_plot = sm.ProbPlot(resid)

pp_qq_plot.ppplot(line = '45')

plt.title('P-P图')

pp_qq_plot.qqplot(line = 'q')

plt.title('Q-Q图')

plt.show()

#残差的正态性检验(非参数法)

standard_resid = (resid-np.mean(resid))/np.std(resid)

g = stats.kstest(standard_resid,'norm')

print(g)

# 总结：由于shapiro正态性检验对样本量的需求是5000以内，而本次数据集样本量有9000多，故选择k-s来完成正态性检验。

# 从k-s检验的p值来看，拒绝了残差服从正态分布的假设，即认为残差并不满足正态性假设这个前提。

# 如果残差不服从正态分布的话，建议对Y变量进行box-cox变换处理。

# 由于fit2模型的残差并没有特别明显的偏态(偏度为0.058，接近于0)，故这里就不对Y进行变换。

#

# import scipy.stats as stats

# #找到box-cox变换的Lambda系数

# lamd = stats.boxcox_normmax(vif.y,method = 'mle')

# #对y进行变换

# vif['trans_y'] = stats.boxcox(vif.y,lamd)

# #建模

# fit3 = sm.formula.ols('y~x1+x2...',data = vif).fit()

# fit3.summary()