#-*- 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()