import pandas as pd

import numpy as np

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 scipy.io

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

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

# 数据读取

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

# data = scipy.io.loadmat('ENCDATA-2hp.mat') # 读取mat文件

# # path = scio.loadmat('FFT-2hp.mat')['FFT-2hp']

# print(data)

# train1=data['train3hp']

# test1=data['test3hp']

# train_y=data['train_y3hp']

# test_y=data['test_y3hp']

# data1=train1[:,1]

# print("train1",train1.shape)

# print("data",data1.shape)

# # sns.pairplot(data)

# # plt.show()

# #y, X = dmatrices( data1, data = train1, return_type= 'dataframe')

# fit2 = sm.formula.ols( data1,data = train1).fit()

# print("fit2",fit2)

# fit2.summary()

# pred2 = fit2.predict()

# print("pred2",pred2)

#

# np.sqrt(mean_squared_error(train1.PE, pred2))

# resid = fit2.resid

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

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

ccpp.describe()

sns.pairplot(ccpp)

plt.show()

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

ccpp.corrwith(ccpp.PE)

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.columnsvif

# print(vif[ "features"])

# 构造PE与AT、V和AP之间的线性模型

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

fit.summary()

print("fit",fit)

# 计算模型的RMSE值

pred = fit.predict()

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

# 离群点检验

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)

ccpp_outliers.head()

print("contat1",contat1)

# 重新建模

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

fit2.summary()

# 计算模型的RMSE值

pred2 = fit2.predict()

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

function(){ //K线图 http://www.kaifx.cn/mt4/kaifx/1770.html

resid = fit2.resid

# 标准化残差与预测值之间的散点图

plt.scatter(fit2.predict(), (fit2.resid-fit2.resid.mean())/fit2.resid.std())

plt.xlabel( '预测值',fontdict={'family' : 'sans-serif', 'size' : 20})

plt.ylabel( '标准化残差',fontdict={'family' : 'sans-serif', 'size' : 20})

# 添加水平参考线

plt.axhline(y = 0, color = 'r', linewidth = 2)

plt.show()