机器学习算法——多元线性回归

数据来源合鲸社区-工业产生的蒸汽量预测

1.1 背景描述

针对当前中央提出碳中和策略，我们提取其中部分关键信息，得知火力发电的基本原理是：燃料在燃烧时加热水生成蒸汽，蒸汽压力推动汽轮机旋转，然后汽轮机带动发电机旋转，产生电能。在这一系列的能量转化中，影响发电效率的核心是锅炉的燃烧效率，即燃料燃烧加热水产生高温高压蒸汽。锅炉的燃烧效率的影响因素很多，包括锅炉的可调参数，如燃烧给量，一二次风，引风，返料风，给水水量；以及锅炉的工况，比如锅炉床温、床压，炉膛温度、压力，过热器的温度等。

经脱敏后的锅炉传感器采集的数据（采集频率是分钟级别），根据锅炉的工况，预测产生的蒸汽量。

环境要求：飞桨 PaddlePaddle 2.2 及以上版本

数据说明 数据分成训练数据（train.txt）和测试数据（test.txt），其中字段”V0”-“V37”，这38个字段是作为特征变量，”target”作为目标变量。选手利用训练数据训练出模型，预测测试数据的目标变量，排名结果依据预测结果的MSE（mean square error）。

模型选择

根据特征变量预测目标变量归类为机器学习中有监督学习回归问题。

1.3 源数据解析

求解回归问题有两种方案： 第一种是用矩阵直接求出w的值，然后做预测就行，但是前提是矩阵不能是奇异矩阵。 第二种是采用梯度下降法求出w和b的局部最优解。

1.3.1 首先采用第一种方案-矩阵法

获取数据

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
df=pd.read_csv("zhengqi_train.txt",sep='\t')
df

image.png
训练集数据由38个特征值和1个目标值构成，首先要做的是将特征值和目标值分开。矩阵运算有个公式： image.png

我们用这个构造一个求w的函数。

def standRegres(df):
    xMat = np.mat(df.iloc[:,:-1].values) # 提取特征
    yMat = np.mat(df.iloc[:,-1].values).T # 提取标签
    xTx = xMat.T*xMat
    if np.linalg.det(xTx) ==0:
        print('这是个奇异矩阵，不能求解逆矩阵') #行列式为0,则该矩阵为奇异矩阵,无法求解逆矩阵
    return
    w = xTx.I * (xMat.T*yMat)
    return w

这里需要提前判断xTx是否是满秩矩阵。若不满秩，则无法对其进行求逆矩阵的操作。定义完函数后，即可测试运行效果。这里需要注意的是，当使用矩阵分解来求解多元线性回归时，必须添加一列全为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)

这是个奇异矩阵，不能求解逆矩阵

由于得出的是奇异矩阵，说明想通过简单的线性代数求解是不行的，只能采用梯度下降的方法来做。

第二种方法：梯度下降法

原始线性方程

向量点乘公式python函数

def predict(x, w, b): 
    p = np.dot(x, w) + b     
    return p 

多元线性回归的损失函数

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

多元梯度下降原理公式

梯度函数

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

梯度下降函数

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

导入我们的训练集得出训练参数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]}")

image.png

通过sklearn库实现

通过sklearn库还是比较方便的，前面的部分都是理解算法过程的，实际编写代码还是通过库比较方便。

from sklearn.linear_model import LinearRegression
reg = LinearRegression()
reg.fit(df1.iloc[:,:-1].values,df1.iloc[:,-1].values)
print(reg.coef_) # 查看系数
reg.intercept_ # 查看截距

image.png

通过判定系数来查看下模型是否准确

####判定系数
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)

image.png

0.889的判定系数还算可以。

最后一步：将测试集放入模型输出预测结果

通过计算好的w和b值将测试集的数据进行拟合得出结果

df2=pd.read_csv("zhengqi_test.txt",sep='\t')
df2
reg.coef_
df2['prdit']=np.dot(df2,reg.coef_)+reg.intercept_
df2

image.png

不管是简单的线性回归或者是多元线性回归对拟合的数据都不是特别的优秀，原因是数据可能是非线性关系，或者可以通过2次或者多次方程拟合更为优秀的模型，但是线性回归是基础。