import pandas as pd
成都创新互联专注为客户提供全方位的互联网综合服务，包含不限于成都网站制作、成都网站设计、延庆网络推广、微信小程序开发、延庆网络营销、延庆企业策划、延庆品牌公关、搜索引擎seo、人物专访、企业宣传片、企业代运营等，从售前售中售后，我们都将竭诚为您服务，您的肯定，是我们最大的嘉奖；成都创新互联为所有大学生创业者提供延庆建站搭建服务，24小时服务热线：13518219792，官方网址：www.cdcxhl.com
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()
