循环发电厂电能输出预测

本次项目是就某联合循环发电厂的数据，运用线性回归模型进行预测电能输出，若文中出现错误的地方，还望指正，谢谢！

1. 数据来源及背景

2. 数据探索分析

3. 相关分析

4. 回归分析

5. 多重共线性

6. 模型应用

1. 数据来源及背景

数据来源: http://archive.ics.uci.edu/ml/machine-learning-databases/00294/

该数据集是收集于联合循环发电厂的9568个数据点, 共包含5个特征: 每小时平均环境变量温度（AT），环境压力（AP），相对湿度（RH），排气真空（V）和净每小时电能输出（PE）, 其中电能输出PE是我们要预测的变量.

2. 数据探索分析

由于我们的数据是excel格式, 而pandas处理excel文件依赖xlrd, 因此我们首先需要安装它.

1. 安装xlrd

pip install xlrd

2. 读取数据

excel文件中共有5个sheet, 我们读取最后一个

import pandas as pd
#第二参数代表要读取的sheet, 0表示第一个, 1表示第二个..., pandas默认读取第一个
df = pd.read_excel(r'D:\Data\CCPP.xlsx', 4)

3. 查看数据前3行/后3行

pd.set_option('display.max_rows', 6)
df

数据维度9568行X5列, 均是数值型数据

4. 查看数据整体信息

df.info()

RangeIndex: 9568 entries, 0 to 9567
Data columns (total 5 columns):
AT    9568 non-null float64
V     9568 non-null float64
AP    9568 non-null float64
RH    9568 non-null float64
PE    9568 non-null float64
dtypes: float64(5)
memory usage: 373.8 KB

无缺失数据, 且均为浮点型.

5. 描述性统计

pd.set_option('display.max_rows', None)
df.describe()

温度(AT): 范围1.81--37.11, 均值19.65, 中位数20.35, 左偏分布

排气真空(V): 范围25.36--81.56, 均值54.31, 中位数52.08, 右偏分布

环境压力(AP): 范围992.89--1033.30, 均值1013.26, 中位数1012.94, 右偏分布

相对湿度(RH): 范围25.56--100.16, 均值73.31, 中位数74.975, 左偏分布

电能输出(PE): 范围420.26--495.76, 均值454.37, 中位数451.55, 右偏分布

通过中位数和均值可以大致看出分布状况, 但是对于具体的偏斜程度, 就需要用偏态系数去衡量

6. 偏态系数

偏态系数反映的是数据的偏斜程度, 若偏态系数大于1或小于-1, 则称为高度偏态分布, 若在0.5~1或-1~-0.5之间, 称为中度偏态分布, 否则, 称为低度偏态分布

for i in df.columns:
    print('{}偏态系数: '.format(i), df[i].skew())

AT偏态系数为:  -0.1363930494749227
V偏态系数为:  0.19852101136676173
AP偏态系数为:  0.2654446935825862
RH偏态系数为:  -0.4318387491833329
PE偏态系数为:  0.3065094354204023

可以看出以上均为低度偏态分布. 我们来看看分布形状的另一个度量--峰态系数

7. 峰态系数

峰态系数反映的是数据的扁平程度, 峰态常常是与标准正态分布对比, 若峰态系数为0, 则说明服从标准正态分布, 若大于0, 则说明比标准正态分布更尖, 称为尖峰分布, 若小于0, 则称为平峰分布.

for i in df.columns:
    print('{}峰态系数为: '.format(i), df[i].kurt())

AT峰态系数为:  -1.0375491923092457
V峰态系数为:  -1.4443366772319615
AP峰态系数为:  0.0942371953033132
RH峰态系数为:  -0.4445263744874084
PE峰态系数为:  -1.0485209686925079

可以看得出: AT, V, AP, PE与正态分布相比较扁平, AP则略尖峰.

 

现在我们已经了解了数据的大致情况, 下面让我们进行下一环节: 相关分析

3. 相关分析

该环节的目的是了解特征与特征之间的依赖程度, 这种依赖程度可以用相关系数来衡量, 而通常在相关系数计算之前, 会通过散点图查看分布状况.

1. 散点图

import seaborn as sns
import matplotlib.pyplot as plt 
#初始化
sns.set()
#绘制分布矩阵
sns.pairplot(df)
#保存图片
plt.savefig('ccpp.png')
plt.show()

可以看得出AT, V, AP, PH分别与PE有不同程度的相关关系, 其中AT, V呈现负相关, AP, PH呈现正相关.

我们可以通过散点图可以大致看出相关关系, 对于关系强度我们则需要用相关系数去度量.

2. 相关系数 

#计算相关系数, 默认为皮尔逊相关系数
correlation = df.corr()
correlation

皮尔逊相关系数反映的是线性相关程度, 其取值介于-1~1之间. 若取0, 说明两变量无关, 若取+1, 说明完全正相关, 若取-1, 说明完全负相关. 

可以将相关系数绘制成热力图

3. 热力图

#绘制热力图, 设置线宽, 最大值, 最小值, 线条白色, 显示数值, 方形
sns.heatmap(correlation, linewidths=0.2, vmax=1, vmin=-1, linecolor='w',
            annot=True,square=True)
plt.savefig('correlation.png')
plt.show()

可以根据热力图的颜色来判断关系强度: 自变量中AT和V与因变量PE为高度相关(0.8~1或-1~-0.8), 自变量AP与因变量PE为中度相关(0.5~0.8或-0.8~-0.5), 自变量RH与因变量PE为低度相关(0.3~0.5或-0.5~-0.3). 另外, 自变量中V和AT高度相关.

接下来我们开始回归分析

4. 回归分析

相关分析是确定是否存在关系以及关系强度, 而回归分析则是确定它们到底是什么样的关系, 建立具体的关系表达式.

由于我们的自变量有4个且各自变量与因变量均有不同程度的线性相关关系, 因此我们假设全部自变量与因变量满足多元线性回归模型:

回归分析的目的就是确定这, , , ,  5个模型参数, 我们这里直接采用sklearn库的线性回归模型(基于最小二乘法)去求解这个5参数.

在求解之前, 将数据集划分为训练集和测试集.

1. 划分数据集

from sklearn.model_selection import train_test_split
X, y = df[['AT', 'V', 'AP', 'RH']], df['PE']
#按照8:2的比例划分为训练数据集和测试数据集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.2)

2. 构建模型

from sklearn.linear_model import LinearRegression
#模型实例化
LR = LinearRegression()
#训练模型
LR.fit(X_train, y_train)
print("截距: ", LR.intercept_)
print("回归系数: ", LR.coef_ )

截距:  445.62291957190826
回归系数:  [-1.97153297 -0.23188593  0.07065982 -0.15715584]

因此, 我们的线性回归模型为: PE=455.62-1.97153297AT-0.23188593V+0.07065982AP -0.15715584RH

3. 模型评估

对于回归模型, 常用的评估方法有: 平均绝对误差(MAE), 均方误差(MSE), 均方根误差(RMSE), 多重判定系数(R²).

1) 平均绝对误差MAE(mean_absolute_error)

其中, 表示实际值, 表示估计值. 平均绝对误差可以有效避免误差相互抵消的问题, 可以更好反映预测误差的实际情况

from sklearn import metrics
#分别对训练集和测试集进行预测
y_train_pred = LR.predict(X_train)
y_test_pred = LR.predict(X_test)
#分别计算训练集和测试集的平均绝对误差
y_train_mae = metrics.mean_absolute_error(y_train, y_train_pred)
y_test_mae = metrics.mean_absolute_error(y_test, y_test_pred)
print('训练集MAE: ', y_train_mae)
print("测试集MAE: ", y_test_mae)

训练集MAE:  3.6280382037797847
测试集MAE:  3.629089312294724

在测试集上预测值与实际值的平均误差比训练集略大.

 

避免误差抵消的另外一种方法就是求平方, 也就是均方误差

 2) 均方误差MSE(mean_square_error)

#分别计算训练集和测试集的均方误差
y_train_mse = metrics.mean_squared_error(y_train, y_train_pred)
y_test_mse = metrics.mean_squared_error(y_test, y_test_pred)
print('训练集MSE: ', y_train_mse)
print("测试集MSE: ", y_test_mse)

训练集MSE:  20.684348049573618
测试集MSE:  21.115806888804244

 同样是测试集比训练集误差大, 由于其数值的单位也是平方级别, 为了更好地描述对其开根号, 也就是均方根误差.

3) 均方根误差RMSE(root_mean_square_error)

from math import sqrt
y_train_rmse = sqrt(y_train_mse)
y_test_rmse = sqrt(y_test_mse)
print('训练集RMSE: ', y_train_rmse)
print("测试集RMSE: ", y_test_rmse)

训练集RMSE:  4.548004842738584
测试集RMSE:  4.5951938902296865

这样就方便去描述了, 即用我们所建立的回归模型来预测每小时净电能输出时, 平均的预测误差为4.6MW, 与训练集相比略大.

 4) 多重判定系数R² (r2_score)

其中, SST则为总误差平方和, 即SST=SSR+SSE, SSR为回归平方和, 反映的是能被回归方程解释的误差; SSE为残差平方和, 反映的是不能被回归方程解释的误差. 也就是说, 多重判定系数反映的是被回归方程解释的误差占总误差的比例, 回归平方和所占比重越大, 说明回归方程拟合程度越好

#分别计算训练集和测试集的多重判定系数
y_train_r2 = metrics.r2_score(y_train, y_train_pred)
y_test_r2 = metrics.r2_score(y_test, y_test_pred)
print('训练集R2: ', y_train_r2)
print("测试集R2: ", y_test_r2)

训练集R2:  0.928600622449398
测试集R2:  0.9289130339796056

即在训练集和测试集上总变差中能被回归方程所解释的比例分别为92.86%, 92.89%, 在测试集上比训练集表现略好.

在sklearn中还有自带的评估器, 不妨来看看.

#直接用训练好的模型去评分
y_train_score = LR.score(X_train, y_train)
y_test_score = LR.score(X_test, y_test)
print('训练集score: ', y_train_score)
print("测试集score: ", y_test_score)

训练集score:  0.928600622449398
测试集score:  0.9289130339796056

可以看得出与多重判定系数的结果是一致的, 说明sklearn自带的评估器采用的是多重判定系数.

由于我们是假设全部自变量与因变量服从线性关系, 那么到底是否真正服从呢? 这就需要进行检验才能得知, 即进入下一步: 模型检验

4. 模型检验

首先, 我们需要检验模型是否服从线性关系,

1) 线性关系检验

由于该检验统计量是以残差平方和(SSE)和回归平方和(SSR)为基础, 故除以它们各自的自由度来构造F检验统计量,

即服从分子自由度为k, 分母自由度为n-k-1的F分布, 其中k为自变量个数, n为样本量. 进行假设检验主要分为5个步骤:

① 提出原假设和备择假设: H₀: , H₁: 不全为0

② 确定显著性水平: 

③ 选择检验统计量: 

④ 建立决策准则:有两种方式一种是拒绝域, 一种是P值.

a. 拒绝域: 即拒绝原假设的区域, 该拒绝域通过临界点 (由于F统计量为二次型, 因此取值非负, 仅有一个拒绝域, 故为α) 来确定, 如果F检验统计量落入拒绝域, 则拒绝原假设从而接受备择假设, 否则不能拒绝原假设.

b. P值: P值表示原假设发生的可能性(概率), 若可能性小于显著性水平, 即小概率事件, 那么就拒绝原假设从而接受备择假设, 否则, 不能拒绝原假设. 

⑤ 下结论: 根据决策准则, 做出是否拒绝原假设的结论.

from pandas import Series
from scipy.stats import f
#将array转为series格式
y_train_pred = Series(y_train_pred, index=y_train.index)
#分别计算训练数据上的SSR和SSE
y_train_ssr = y_train_pred.apply(lambda x: (x - y_train.mean())**2).sum()
y_train_sse = y_train.sub(y_train_pred).apply(lambda x: x**2).sum()
#dn是SSR的自由度(自变量个数), df则是SSE的自由度
dn, df = 4, y_train.shape[0]-dn-1
#计算F值
y_train_f = (y_train_ssr/dn) / (y_train_sse/df)
#计算p值
p = f.sf(y_train_f, dn, df)
#计算0.05显著性水平下临界值
cr_value = f.isf(0.05, dn, df)
print('训练数据集F值: ', y_train_f)
print('0.05显著性水平下临界值: ', cr_value)
print('训练数据集P值: %.20f'% p)

训练数据集F值:  24870.196368594174
0.05显著性水平下临界值:  2.3730935370191646
训练数据集P值: 0.00000000000000000000

F值大于临界值, 拒绝原假设H₀ ( P值小于显著性水平, 同样可以拒绝原假设), 表明全部自变量与因变量的线性关系是显著的, 但这并不意味着每个自变量都和因变量的线性关系显著, 因此还需要进行回归系数(自变量)检验.

2) 回归系数检验

对每个回归系数构造t检验统计量:

其中, 表示第i个回归系数的估计值(因为是根据最小二乘法估计的, 而不是真实值), 表示第i个回归系数的标准误差. 则:

其中, MSE表示无偏残差平方和, 即残差平方和SSE除以它的自由度(n-k-1), C_ii表示(XX^T)^-1逆矩阵i行第i列的元素. t检验过程同样是五步骤:

① 提出原假设和备择假设: H₀: , H₁: 不为0

② 确定显著性水平: 

③ 选择检验统计量: 

④ 建立决策准则: 同样有两种方式, 不过t检验的拒绝域有两个.

⑤ 下结论: 根据决策准则做出是否拒绝原假设的结论

from scipy.stats import t
def get_tvalue(sse, df, matr, beta, i):
    '''计算t值'''
    mse = sse / df
    sbeta = sqrt(matr[i+1, i+1]* mse)
    t = beta / sbeta
    return t

limit = t.isf(0.025, df)
print('0.05显著性水平下的临界值: ', limit)

X_train['B'] = 1
X_train = X_train.reindex(columns=['B', 'AT', 'V', 'AP', 'RH'])
#转成矩阵
xm = np.mat(X_train)
#计算(X'X)的逆矩阵
xmi = np.dot(xm.T, xm).I
index, betas = range(4), LR.coef
for i, beta in zip(index, betas):
    tvalue = get_tvalue(y_train_sse, df, xmi, beta, i)
    pvalue = t.sf(abs(tvalue), df)*2
    print('beta{0}的t值: '.format(i+1), tvalue)
    print('beta{0}的p值: '.format(i+1), pvalue)

0.05显著性水平下的临界值:  2.085963447265837
beta1的t值:  -5.898890146924707
beta1的p值:  9.049074948455422e-06
beta2的t值:  -1.4562620735727096
beta2的p值:  0.16084038317451793
beta3的t值:  0.34176894237398453
beta3的p值:  0.7360897563076247
beta4的t值:  -1.734158669159414
beta4的p值:  0.09827823435663764

从以上结果中可以看到只有beta1(自变量AT)的t值的绝对值大于临界值, 即只有beta1拒绝原假设, 也就是说4个自变量中只有beta1影响是显著的, 那么为什么其他自变量不呢? 造成这种的原因可能是自变量之间高度相关, 这种问题称为多重共线性.

5. 多重共线性

引用百度百科的定义

多重共线性是指线性回归模型中的解释变量之间由于存在精确相关关系或高度相关关系而使模型估计失真或难以估计准确。

根据定义可知, 多重共线性的判别方法是相关系数矩阵, 查看之前的相关系数矩阵发现在AT和V这两个解释变量(自变量)之间存在高度相关关系: 数值为0.84. 因此可以认为变量间存在多重共线性.

对于多重共线性的处理方法有: 删除共线变量, 转换模型的形式, 逐步回归法, 岭回归法, 主成分分析法. 在这里我采用前两种方法.

1. 删除共线变量: 通常是删除高度相关的变量, 保留重要的变量. 重要的变量是指通过t检验的变量, 于是我们删除V这个解释变量.

X1, y1 = df[['AT', 'AP', 'RH']], df['PE']
X1_train, X1_test, y1_train, y1_test = train_test_split(X1, y1, test_size = 0.2)
#训练模型, 预测以及计算均方根误差和多重判定系数
LR1 = LinearRegression()
LR1.fit(X1_train, y1_train)
inter, co = LR1.intercept_, LR1.coef_
y1_train_pred = LR1.predict(X1_train)
y1_train_rmse = sqrt(metrics.mean_squared_error(y1_train, y1_train_pred))
y1_train_score = LR1.score(X1_train, y1_train)
print("回归模型: PE={0}+{1}AT+{2}AP+{3}RH".format(inter,co[0], co[1], co[2]))
print('训练集RMSE: ', y1_train_rmse)
print('训练集拟合优度: ', y1_train_score)

#计算F检验中0.05显著性水平下的P值
y1_train_pred = Series(y1_train_pred, index=y1_train.index)
y1_train_ssr = y1_train_pred.apply(lambda x: (x - y1_train.mean())**2).sum()
y1_train_sse = y1_train.sub(y1_train_pred).apply(lambda x: x**2).sum()
dn1, df1 = 3, y1_train.shape[0]-3-1
y1_train_f = (y1_train_ssr/dn1) / (y1_train_sse/df1)
y1_p = f.sf(y1_train_f, dn1, df1)
# cr_value = f.isf(0.05, dn, df)
print('F检验 0.05显著性水平下训练数据的P值: %.20f'% y1_p)

#计算t检验在0.05显著性水平下的P值
def get_t1value(sse, df, matr, beta, i):
    mse = sse / df
    sbeta = sqrt(matr[i+1, i+1]* mse)
    t = beta / sbeta
    return t
X1_train['B'] = 1
X1_train = X1_train.reindex(columns=['B', 'AT', 'AP', 'RH'])
xm1 = np.mat(X1_train)
xmi1 = np.dot(xm1.T, xm1).I
index, betas = range(3), LR1.coef_
for i, beta in zip(index, betas):
    tvalue = get_t1value(y1_train_sse, df1, xmi1, beta, i)
    pvalue = t.sf(abs(tvalue), df1)*2
    print('t检验 0.05显著性水平下beta{0}的p值: '.format(i+1), pvalue)

View Code

回归模型: PE=485.39627734198206+-2.374833922164532AT+0.030214741968456443AP+-0.20456036084843196RH
训练集RMSE:  4.839764102773314
训练集拟合优度:  0.9195738379930241
F检验 0.05显著性水平下训练数据的P值: 0.00000000000000000000
t检验 0.05显著性水平下beta1的p值:  0.0
t检验 0.05显著性水平下beta2的p值:  0.006605096641032901
t检验 0.05显著性水平下beta3的p值:  0.0

均方根误差比之前增大0.3, 拟合优度相比之前减少0.009, F检验通过, 且回归系数均通过, 可见删除共线变量还是可以有效解决多重共线性问题, 但是误差相比之前增大0.2, 因此这种方法还是有一定局限性. 换另外一种方法试试.

2. 转换模型形式: 将数据进行对数转换.

df_log = np.log(df)

回归模型: PE=2.8675586960921957+-0.05125235913516803AT+-0.05608345239278856V+0.5268863140584524AP+-0.0059335463238909605RH
训练集RMSE:  0.010378041999902672
训练集拟合优度:  0.9229456762426529
F检验 0.05显著性水平下训练数据的P值: 0.00000000000000000000
t检验 0.05显著性水平下beta1的p值:  0.0
t检验 0.05显著性水平下beta2的p值:  0.0
t检验 0.05显著性水平下beta3的p值:  8.936059858590797e-108
t检验 0.05显著性水平下beta4的p值:  1.7118836102396494e-20

同样也都通过检验, 而且经过对数转换后的均方根误差更小. 对比可以发现进行对数转换比删除变量的拟合优度更好些. 我们接下来便是用转换后的模型去预测了.

 6. 模型应用

import matplotlib.pyplot as plt
#用训练好的模型去预测
y2_test_pred = LR2.predict(X2_test)
y2_test_rmse = sqrt(metrics.mean_squared_error(y2_test, y2_test_pred))
# y1_test_rmse = sqrt(metrics.mean_squared_error(y1_test, y1_test_pred))
y2_test_score = LR2.score(X2_test, y2_test)
print('测试集RMSE: ', y2_test_rmse)
print('测试集拟合优度: ', y2_test_score)
#绘制曲线
plt.figure(figsize=(12,8))
plt.plot(range(len(y2_test)), y2_test, 'g', label='test data')
plt.plot(range(len(y2_test_pred)), y2_test_pred, 'r', label='predict data', linewidth=2, alpha=0.8)
plt.legend()
plt.savefig('tp.png')
plt.show()

测试集RMSE:  0.01038585831215862
测试集拟合优度:  0.9223780289364194

可以看到在测试集上表现还是相当不错的.

以上便是对线性回归模型学习过程的总结, 若存在问题, 希望指出, 谢谢!

参考

网易云课堂《吴恩达机器学习》

https://blog.csdn.net/hhsh49/article/details/53169121

https://www.cnblogs.com/pinard/p/6016029.html