线性回归

线性回归(linear regression)是一种建立线性预测函数的模型，通过已有的数据集估计模型的参数，建立自变量和因变量之间的关系。

线性模型
从机器学习的角度看，自变量为数据集X,因变量为y,假设的线性函数模型为
$f(x;w,b)=w^Tx+b$
其中w,b都是可学习的参数，一般为了简化可以将b加入到w的参数向量中，x中添加一个1的特征，所以模型可以简化为
$f(x;w)=w^Tx$

2.参数学习
线性模型中参数学习方法为最小二乘法，线性回归的模型输出值和真实值都为实数，可以使用平方损失函数来衡量真实值和预测值的差异。首先将其转换为数学表达式，数据集D={(x1,y1),(x2,y2)......(xn,yn)},如果每一个输入数据x有d维特征，一般情况x的输入向量默认为列向量。数据集D的输入矩阵表示如下
$X=(x_{1},x_{2},...,x_{n})=\left[ \begin{matrix} x_{1}^T \\ x_{2}^T \\ . \\ .\\ .\\x_{n}^T \end{matrix} \right]= \left[ \begin{matrix} x_{11}^T&x_{12}^T& .& .& .&x_{1d}^T \\ x_{21}^T& x_{22}^T&.&.&. &x_{2d}^T\\ . \\ .\\ .\\x_{n1}^T &x_{n2}^T&.&.&.&x_{nd}^T\end{matrix} \right]$
输出的y表示
$y=\left[ \begin{matrix} y_{1} \\ y_{2}\\ . \\ .\\ .\\y_{n}\end{matrix} \right]$
最小二乘的损失函数如下：
$L(W)=\frac{1}{2m}\sum_{i=1}^{n}(y^n-w^Tx^n)\\ =\frac{1}{2m}\left\|y-XW\right\|_2$
对w求偏导
$\frac {\partial L(W)}{\partial w}=\frac{1}{2m}\frac{\partial \left\|y-XW\right\|_2}{\partial w}\\ =-X^T(y-XW)$
令 $\frac {\partial L(W)}{\partial w}=0$ 得到的最优参数 $W^*$ 为
$W^*={(X^TX)}^{-1}X^Ty\\$

3.从概率角度看最小二乘法
假设y是一个随机变量服从均值为 $W^TX$ ,方差为 ${\delta}^2$ 的高斯分布。
$p(y|x;w,\delta)=N(y;W^TX,{\delta}^2)\\ =\frac{1}{\delta\sqrt{2\pi}}exp(-\frac{y-w^Tx}{2\delta^2})$
参数w在数据集D上的极大似然估计(Likelihood)为
$p(y|X;w,\delta)=\prod_{i=1}^{n}p(y^n|x^n;w,\delta)$
为了方便计算，对似然函数取对数似然函数
$logp(y|X;w,\delta)=log\prod_{i=1}^{n}p(y^i|x^i;w,\delta)\\ =\sum_{i=1}^{n}logp(y^i|x^i;w,\delta)\\ =\sum_{i=1}^{n}log\frac{1}{\delta\sqrt{2\pi}}exp(-\frac{y^i-w^Tx^i}{2\delta^2})\\ =\sum_{i=1}^{n}log\frac{1}{\delta\sqrt{2\pi}}+\sum_{i=1}^{n}logexp(-\frac{y^i-w^Tx^i}{2\delta^2})\\ =\sum_{i=1}^{n}(log\frac{1}{\delta\sqrt{2\pi}}-\frac{y^i-w^Tx^i}{2\delta^2})$
令 $\frac{\partial logp(y|X;w,\delta)}{\partial w}=0$ 得到
$w^{ML}=(X^TX)^{-1}X^Ty\\$
可以看出极大似然解和最小二乘的解相同

4.最小二乘法加入正则化(岭回归）
从最小二乘法的解 $w*=(X^TX)^{-1}X^Ty$ 的公式中看出如果X^TX不可逆，就无法求出w参数，为了解决这个问题提出了岭回归可以加入一个常数 $\lambda$ 的单位矩阵，使其行列式值不为0。参数w*为
$w*=(X^TX+\lambda I)^{-1}X^Ty\\$
岭回归可以看作是加入了正则化的最小二乘估计
$L(W)=\frac{1}{2m}\left\|y-XW\right\|_2+\frac{1}{2}\lambda{\left\|w\right\|_2}^2\\$

5.从概率角度看加入正则化（最大后验估计）
从贝叶斯的角度看参数w不是常向量，而是一个随机向量服从一个先验分布 $p(w;\delta_{0})$ ,一般令w服从0均值的高斯分布
$p(w;\delta_{0})=N(0,{\delta_{0}}^2)\\$
根据贝叶斯公式，w的后验分布
$p(w|X,y;\delta_{0},\delta)\propto p(y|X;w,\delta)p(w,\delta_{0})\\$
最大后验估计（Maximum A Posteriori Estimation，MAP）是指最优参数为后验分布
$p(w|X,y;\delta_{0},\delta)$ 中概率密度最高的参数w
$w^{MAP}=argmaxp(y|X;w,\delta)p(w,\delta_{0})\\ =argmaxlogp(y|X;w,\delta)p(w,\delta_{0})\\ =argmaxlog(\frac{1}{\delta\sqrt{2\pi}}\frac{1}{\delta_{0}\sqrt{2\pi}}-\frac{{(y-w^Tx)}^2}{2\delta^2}-\frac{{\left\|w\right\|_2}^2}{2\delta_{0}^2}) =argmax(-\frac{{(y-w^Tx)}^2}{2\delta^2}-\frac{{\left\|w\right\|_2}^2}{2\delta_{0}^2})$
可以看出最大后验概率估计等价于加入正则项的最小二乘法

6.线性回归python 代码
用sklearn实现线性回归特别简单

import matplotlib.pyplot as plt
import pandas as pd
import numpy as np


data = pd.read_csv('ccpp.csv')
X = data[['AT', 'V', 'AP', 'RH']]
y = data[['PE']]

from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=1)

from sklearn.linear_model import LinearRegression

linreg = LinearRegression()

linreg.fit(X_train, y_train)

print(linreg.intercept_)
print(linreg.coef_)

y_pred = linreg.predict(X_test)
from sklearn import metrics

print("MSE:", metrics.mean_squared_error(y_test, y_pred))
print("RMSE", np.sqrt(metrics.mean_squared_error(y_test, y_pred)))

fig, ax = plt.subplots()

ax.scatter(y_test, y_pred)
ax.plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()], 'k--', lw=4)
ax.set_xlabel('Measured')
ax.set_ylabel('Predicted')
plt.show()

运行后的结果

### w的参数值
[447.06297099]
[[-1.97376045 -0.23229086  0.0693515  -0.15806957]]
MSE: 20.08040120207389
RMSE 4.481116066570235

linear.png

sklearn中的线性回归是使用最小二乘实现的，可以手写一个python 最小二乘练习一下，看看有没有错

import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn import metrics


def read_csv():
    data = pd.read_csv('ccpp.csv')
    X = data[['AT', 'V', 'AP', 'RH']]
    y = data[['PE']]
    return X, y


def lse(X_train, y_train):
    X = X_train.values
    X = np.insert(X, 0, 1, axis=1)
    y = y_train.values
    XTX = X.T.dot(X)
    if np.linalg.det(XTX) == 0.0:
        print("This matrix is sigular, cannot do inverse")  # 行列式的值为0，无逆矩阵
        return
    w = np.linalg.inv(XTX).dot(X.T.dot(y))
    return w


def predict(X_test, w):
    X = X_test.values
    X = np.insert(X, 0, 1, axis=1)
    return X.dot(w)


if __name__ == '__main__':
    X, y = read_csv()
    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=1)
    w = lse(X_train, y_train)
    print(w)
    y_pred = predict(X_test, w)
    print("MSE:", metrics.mean_squared_error(y_test, y_pred))
    print("RMSE", np.sqrt(metrics.mean_squared_error(y_test, y_pred)))
    fig, ax = plt.subplots()
    ax.scatter(y_test, y_pred)
    ax.plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()], 'k--', lw=4)
    ax.set_xlabel('Measured')
    ax.set_ylabel('Predicted')
    plt.show()

运行后的结果是和sklearn中一样的

#w的参数
[[ 4.47062971e+02 -1.97376045e+00 -2.32290859e-01  6.93514994e-02
  -1.58069568e-01]]
MSE: 20.0804012023701
RMSE 4.481116066603286

如果用最小二乘无法直接求出解析解的，也可以使用梯度下降算法求最优解

import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn import metrics


def read_csv():
    data = pd.read_csv('ccpp.csv')
    X = data[['AT', 'V', 'AP', 'RH']]
    y = data[['PE']]
    return X, y


def compute_gradient(X_train, y_train, w):
    n = X_train.shape[0]
    X = X_train.values
    y = y_train.values
    w_gradient = X.T.dot(X.dot(w) - y) / n
    return w_gradient


def compute_error(w, X_train, y_train):
    n = y_train.shape[0]
    X = X_train.values
    y = y_train.values
    error = np.dot((X.dot(w) - y).T, (X.dot(w) - y)) / (2 * n)
    return error[0][0]


def optimizer(X_train, y_train, learning_rate, num_iter):
    m = X_train.shape[1]
    w = np.zeros((m, 1))
    errors = []
    for i in range(num_iter):
        w_gradient = compute_gradient(X_train, y_train, w)
        w = w - learning_rate * w_gradient
        if i % 500 == 0:
            error = compute_error(w, X_train, y_train)
            # print(i, error)
            errors.append(error)
    return w, errors


def predict(X_test, w):
    X = X_test.values
    return X.dot(w)


if __name__ == '__main__':
    X, y = read_csv()
    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=1)
    X_train.insert(0, 'one', 1)
    X_test.insert(0, 'one', 1)
    w, errors = optimizer(X_train, y_train, learning_rate=0.000001, num_iter=10000)
    y_pred = predict(X_test, w)
    print("MSE:", metrics.mean_squared_error(y_test, y_pred))
    print("RMSE", np.sqrt(metrics.mean_squared_error(y_test, y_pred)))
    fig, (ax1, ax2) = plt.subplots(2, 1, sharey=False)
    ax1.scatter(y_test, y_pred)
    ax1.plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()], 'k--', lw=4)
    ax1.set_xlabel('Measured')
    ax1.set_ylabel('Predicted')
    ax2.plot(np.arange(0, 10000, 500), errors, 'b')
    ax2.set_xlabel('num_iter')
    ax2.set_ylabel('error')
    plt.show()

linear1.png

代码,数据集参考博文
https://www.cnblogs.com/pinard/p/6016029.html

线性回归

推荐阅读更多精彩内容

友情链接更多精彩内容