數(shù)據(jù)來源合鯨社區(qū)-工業(yè)產(chǎn)生的蒸汽量預(yù)測
1.1 背景描述
針對當(dāng)前中央提出碳中和策略，我們提取其中部分關(guān)鍵信息，得知火力發(fā)電的基本原理是：燃料在燃燒時(shí)加熱水生成蒸汽，蒸汽壓力推動汽輪機(jī)旋轉(zhuǎn)，然后汽輪機(jī)帶動發(fā)電機(jī)旋轉(zhuǎn)，產(chǎn)生電能。在這一系列的能量轉(zhuǎn)化中，影響發(fā)電效率的核心是鍋爐的燃燒效率，即燃料燃燒加熱水產(chǎn)生高溫高壓蒸汽。鍋爐的燃燒效率的影響因素很多，包括鍋爐的可調(diào)參數(shù)，如燃燒給量，一二次風(fēng)，引風(fēng)，返料風(fēng)，給水水量；以及鍋爐的工況，比如鍋爐床溫、床壓，爐膛溫度、壓力，過熱器的溫度等。
經(jīng)脫敏后的鍋爐傳感器采集的數(shù)據(jù)（采集頻率是分鐘級別），根據(jù)鍋爐的工況，預(yù)測產(chǎn)生的蒸汽量。
環(huán)境要求：飛槳 PaddlePaddle 2.2 及以上版本
數(shù)據(jù)說明 數(shù)據(jù)分成訓(xùn)練數(shù)據(jù)（train.txt）和測試數(shù)據(jù)（test.txt），其中字段”V0”-“V37”，這38個(gè)字段是作為特征變量，”target”作為目標(biāo)變量。選手利用訓(xùn)練數(shù)據(jù)訓(xùn)練出模型，預(yù)測測試數(shù)據(jù)的目標(biāo)變量，排名結(jié)果依據(jù)預(yù)測結(jié)果的MSE（mean square error）。
模型選擇
根據(jù)特征變量預(yù)測目標(biāo)變量歸類為機(jī)器學(xué)習(xí)中有監(jiān)督學(xué)習(xí)回歸問題。
1.3 源數(shù)據(jù)解析
求解回歸問題有兩種方案： 第一種是用矩陣直接求出w的值，然后做預(yù)測就行，但是前提是矩陣不能是奇異矩陣。 第二種是采用梯度下降法求出w和b的局部最優(yōu)解。
1.3.1 首先采用第一種方案-矩陣法
獲取數(shù)據(jù)
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
df=pd.read_csv("zhengqi_train.txt",sep='\t')
df
訓(xùn)練集數(shù)據(jù)由38個(gè)特征值和1個(gè)目標(biāo)值構(gòu)成，首先要做的是將特征值和目標(biāo)值分開。矩陣運(yùn)算有個(gè)公式：
我們用這個(gè)構(gòu)造一個(gè)求w的函數(shù)。
def standRegres(df):
xMat = np.mat(df.iloc[:,:-1].values) # 提取特征
yMat = np.mat(df.iloc[:,-1].values).T # 提取標(biāo)簽
xTx = xMat.T*xMat
if np.linalg.det(xTx) ==0:
print('這是個(gè)奇異矩陣，不能求解逆矩陣') #行列式為0,則該矩陣為奇異矩陣,無法求解逆矩陣
return
w = xTx.I * (xMat.T*yMat)
return w
這里需要提前判斷xTx是否是滿秩矩陣。若不滿秩，則無法對其進(jìn)行求逆矩陣的操作。定義完函數(shù)后，即可測試運(yùn)行效果。這里需要注意的是，當(dāng)使用矩陣分解來求解多元線性回歸時(shí)，必須添加一列全為1的列，用于表征線性方程截距。
測試一下
ex = pd.DataFrame(np.ones([df.shape[0],1]))
ex
X=df.iloc[:,:-1]
Y=df.iloc[:,-1]
df = pd.concat([ex,X,Y],axis=1)
ws = standRegres(df)
這是個(gè)奇異矩陣，不能求解逆矩陣
由于得出的是奇異矩陣，說明想通過簡單的線性代數(shù)求解是不行的，只能采用梯度下降的方法來做。
第二種方法：梯度下降法
原始線性方程
向量點(diǎn)乘公式python函數(shù)
def predict(x, w, b):
p = np.dot(x, w) + b
return p
多元線性回歸的損失函數(shù)
def compute_cost(X, y, w, b):
m = X.shape[0]
cost = 0.0
for i in range(m):
f_wb_i = np.dot(X[i], w) + b
cost = cost + (f_wb_i - y[i])**2
cost = cost / (2 * m)
return cost
多元梯度下降原理公式
梯度函數(shù)
def compute_gradient(X, y, w, b):
"""
Computes the gradient for linear regression
Args:
X (ndarray (m,n)): Data, m examples with n features
y (ndarray (m,)) : target values
w (ndarray (n,)) : model parameters
b (scalar) : model parameter
Returns:
dj_dw (ndarray (n,)): The gradient of the cost w.r.t. the parameters w.
dj_db (scalar): The gradient of the cost w.r.t. the parameter b.
"""
m,n = X.shape #(number of examples, number of features)
dj_dw = np.zeros((n,))
dj_db = 0.
for i in range(m):
err = (np.dot(X[i], w) + b) - y[i]
for j in range(n):
dj_dw[j] = dj_dw[j] + err * X[i, j]
dj_db = dj_db + err
dj_dw = dj_dw / m
dj_db = dj_db / m
return dj_db, dj_dw
梯度下降函數(shù)
def gradient_descent(X, y, w_in, b_in, cost_function, gradient_function, alpha, num_iters):
"""
Performs batch gradient descent to learn theta. Updates theta by taking
num_iters gradient steps with learning rate alpha
Args:
X (ndarray (m,n)) : Data, m examples with n features
y (ndarray (m,)) : target values
w_in (ndarray (n,)) : initial model parameters
b_in (scalar) : initial model parameter
cost_function : function to compute cost
gradient_function : function to compute the gradient
alpha (float) : Learning rate
num_iters (int) : number of iterations to run gradient descent
Returns:
w (ndarray (n,)) : Updated values of parameters
b (scalar) : Updated value of parameter
"""
# An array to store cost J and w's at each iteration primarily for graphing later
J_history = []
w = copy.deepcopy(w_in) #avoid modifying global w within function
b = b_in
for i in range(num_iters):
# Calculate the gradient and update the parameters
dj_db,dj_dw = gradient_function(X, y, w, b) ##None
# Update Parameters using w, b, alpha and gradient
w = w - alpha * dj_dw ##None
b = b - alpha * dj_db ##None
# Save cost J at each iteration
if i<100000: # prevent resource exhaustion
J_history.append( cost_function(X, y, w, b))
# Print cost every at intervals 10 times or as many iterations if < 10
if i% math.ceil(num_iters / 10) == 0:
print(f"Iteration {i:4d}: Cost {J_history[-1]:8.2f} ")
return w, b, J_history #return final w,b and J history for graphing
導(dǎo)入我們的訓(xùn)練集得出訓(xùn)練參數(shù)w和b
import copy, math
df1=pd.read_csv("zhengqi_train.txt",sep='\t')
w_init=np.array(Y)
X_train=np.array(df1.iloc[:,:-1]).T
y_train=np.array(df1.iloc[:,-1])
df1
# initialize parameters
initial_w = np.zeros_like(w_init)
initial_b = 0.
# some gradient descent settings
iterations = 1000
alpha = 5.0e-7
# run gradient descent
w_final, b_final, J_hist = gradient_descent(X_train, y_train, initial_w, initial_b,
compute_cost, compute_gradient,
alpha, iterations)
print(f"b,w found by gradient descent: {b_final:0.2f},{w_final} ")
m,_ = X_train.shape
for i in range(m):
print(f"prediction: {np.dot(X_train[i], w_final) + b_final:0.2f}, target value: {y_train[i]}")
通過sklearn庫實(shí)現(xiàn)
通過sklearn庫還是比較方便的，前面的部分都是理解算法過程的，實(shí)際編寫代碼還是通過庫比較方便。
from sklearn.linear_model import LinearRegression
reg = LinearRegression()
reg.fit(df1.iloc[:,:-1].values,df1.iloc[:,-1].values)
print(reg.coef_) # 查看系數(shù)
reg.intercept_ # 查看截距
通過判定系數(shù)來查看下模型是否準(zhǔn)確
####判定系數(shù)
from sklearn.metrics import mean_squared_error,r2_score
yhat = reg.predict(df1.iloc[:,:-1])
yhat
y = df1.iloc[:,-1].values
mean_squared_error(y,yhat)
r2_score(y,yhat)
0.889的判定系數(shù)還算可以。
最后一步：將測試集放入模型輸出預(yù)測結(jié)果
通過計(jì)算好的w和b值將測試集的數(shù)據(jù)進(jìn)行擬合得出結(jié)果
df2=pd.read_csv("zhengqi_test.txt",sep='\t')
df2
reg.coef_
df2['prdit']=np.dot(df2,reg.coef_)+reg.intercept_
df2
不管是簡單的線性回歸或者是多元線性回歸對擬合的數(shù)據(jù)都不是特別的優(yōu)秀，原因是數(shù)據(jù)可能是非線性關(guān)系，或者可以通過2次或者多次方程擬合更為優(yōu)秀的模型，但是線性回歸是基礎(chǔ)。
