#-*- 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
 
#數(shù)據(jù)獲取
ccpp = pd.read_excel('CCPP.xlsx')
ccpp.describe()
#繪制各變量之間的散點(diǎn)圖
sns.pairplot(ccpp)
plt.show()
#發(fā)電量（PE）與自變量之間的相關(guān)系數(shù)
a = ccpp.corrwith(ccpp.PE)
print(a)
#將因變量PE，自變量AT,V,AP和截距項(xiàng)（值為1的1維數(shù)值）以數(shù)據(jù)框的形式組合起來
y,x = dmatrices('PE~AT+V+AP',data = ccpp,return_type = 'dataframe')
#構(gòu)造空的數(shù)據(jù)框
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)
 
#構(gòu)建PE與AT,V和AP之間的線性模型
fit = sm.formula.ols('PE~AT+V+AP',data=ccpp).fit()
b = fit.summary()
# print(b)
#計(jì)算模型的RMSE值
pred = fit.predict()
c = np.sqrt(mean_squared_error(ccpp.PE,pred))
print(c)
#離群點(diǎn)檢驗(yàn)
outliers = fit.get_influence()
#高杠桿值點(diǎn)（帽子矩陣）
leverage = outliers.hat_matrix_diag
#dffits值
dffits = outliers.dffits[0]
#學(xué)生化殘差
resid_stu = outliers.resid_studentized_external
#cook距離
cook = outliers.cooks_distance[0]
#covratio值
covratio = outliers.cov_ratio
#將上面的幾種異常值檢驗(yàn)統(tǒng)計(jì)量與原始數(shù)據(jù)集合并
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)
 
#計(jì)算異常值數(shù)量的比例
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)
#
#殘差的正態(tài)性檢驗(yàn)（直方圖法）
resid = fit2.resid
#中文和負(fù)號的正常顯示
# 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')
#設(shè)置坐標(biāo)軸標(biāo)簽和標(biāo)題
plt.title('殘差直方圖')
plt.ylabel('密度值')
#生成正態(tài)曲線的數(shù)據(jù)
x1 = np.linspace(resid.min(),resid.max(),1000)
normal = mlab.normpdf(x1,resid.mean(),resid.std())
#繪制正態(tài)分布曲線
plt.plot(x1,normal,'r-',linewidth = 2,label = '正態(tài)分布曲線')
#生成核密度曲線的數(shù)據(jù)
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()
#生成的正態(tài)曲線的數(shù)據(jù)
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()
#殘差的正態(tài)性檢驗(yàn)（非參數(shù)法）
standard_resid = (resid-np.mean(resid))/np.std(resid)
g = stats.kstest(standard_resid,'norm')
print(g)
# 總結(jié)：由于shapiro正態(tài)性檢驗(yàn)對樣本量的需求是5000以內(nèi)，而本次數(shù)據(jù)集樣本量有9000多，故選擇k-s來完成正態(tài)性檢驗(yàn)。
# 從k-s檢驗(yàn)的p值來看，拒絕了殘差服從正態(tài)分布的假設(shè)，即認(rèn)為殘差并不滿足正態(tài)性假設(shè)這個(gè)前提。
# 如果殘差不服從正態(tài)分布的話，建議對Y變量進(jìn)行box-cox變換處理。
# 由于fit2模型的殘差并沒有特別明顯的偏態(tài)（偏度為0.058，接近于0），故這里就不對Y進(jìn)行變換。
 
# 
# import scipy.stats as stats
# #找到box-cox變換的Lambda系數(shù)
# lamd = stats.boxcox_normmax(vif.y,method = 'mle')
# #對y進(jìn)行變換
# vif['trans_y'] = stats.boxcox(vif.y,lamd)
# #建模
# fit3 = sm.formula.ols('y~x1+x2...',data = vif).fit()
# fit3.summary()