多元线性回归分析--推导过程及python、numpy和python实现

多元线性回归分析

什么是线性回归？

线性回归，如上图所示（这里用二维的例子比较好理解），我们知道许多的 $(x_1,y_1),(x_2,y_2),...,(x_n,y_n)$，即图中红色的点，通过某种方法，得到图中蓝色的线（$y=w\times x+b$），即求 $w，b$的值；然后可以使得未知数据 $x_{new}$ 输入方程，使得结果 $y_{new}$ 尽可能的接近真实值。

具体的关于线性回归的定义，这里是引用西瓜书上定义

定义：

​ 给定由 $d$ 个属性描述的示例 $x=(x_1;x_1;...;x_d)$，其中 $x_i$ 是 $\text{x}$ 在第 $i$ 个属性上的取值，线性模型（linear model）试图学得一个通过属性的组合来进行预测的函数，即
$$f(x)=w_1{x}_1+w_2{x}_2+...+w_d{x}_d+b$$
一般用向量形式写成
$$f(x)=w^{T}x+b$$
其中 $w=(w_1;w_2;...;w_1;)$。$w$ 和 $b$学习到以后，模型就得以确定

通常回归分为线性回归分析和非线性回归分析，判读二者最好的方法是看每一个变量的指数，如都是1次的就是线性，否则是非线性。

基于矩阵求解（解析解法或者公式法）

推导过程

假如我们的多元线性方程有n项，目标方程为
$$f\left(x_i\right)=w_1{x}_1+w_2{x}_2+\cdots +w_n{x}_n+b$$
有对应的关于数据D（$(x_1,x_2,...,x_n,y)$）,需要根据输入的数据来求解（$w_i,b$）。这里对于最后的偏置项$b$，可以认为是$w_0{x}_0$，其中$x_0$是一个恒为1的数，求解$w_0$相当于求解$b$，即
$$f\left(x_i\right)=w_0{x}_0+w_1{x}_1+w_2{x}_2+\cdots +w_n{x}_n$$
这样写成矩阵的形式为
$$f\left(x\right)=w^{T}x$$
其中$x$是一个列向量，$w$是一个列向量

在数据D中有多组数据，就可以写出矩阵的形式

$$w^T\left[\begin{array}{c}{x}^1,x^2,\dots ,x^n\end{array}\right]=\left[\begin{array}{c}{y}^1,y^2,\dots ,y^n\end{array}\right]$$

如何确定$w$，就是在于如何可以使得$y^i$(目标值)和${\hat{y}}^i$(预测值)直接最小，这里使用均方差来使得误差最小化，即定义下面的目标方程
$$\begin{array}{lc}{w}^{\ast }& =\underset{w}{\mathrm{argmin}}\sum _{i=1}^m(y^i-{\hat{y}}^i{)}^2\\ & =\underset{w}{\mathrm{argmin}}\sum _{i=1}^m(y^i-w^T{x}^i{)}^2\end{array}$$
其中$w^{\ast }$就是我们要找的最优的参数，$m$表示有m组数据。根据上式可以知道，目标函数是一个关于$w$的凸函数，对于凸函数而言，一阶导数等于零对应的$w$就是目标函数的最优值；为了求导方便，我们再把式子变换一下，式子中$x^i$是一个列向量，如果有多组数据的话就可以把他写成一个矩阵$X=(x^1,x^2,...,x^m)$，相应的真实值$y=(y^1,y^2,...,y^n{)}^T$是一个列向量，然后目标函数就可以转换为
$$E=(y-X^{T}w{)}^T(y-X^{T}w)$$
对$w$求导
$$\frac{\partial E}{\partial w}=2X(X^{T}w-y)$$
令$\frac{\partial E}{\partial w}=0$，就可以求的
$$w=(X{X}^T{)}^{-1}Xy$$

代码实现

import numpy as np

def getParams(x, y):
    n = len(y)
    x_0 = np.ones(n)
    x = np.matrix(np.vstack([x_0, x]).T)
    y = np.matrix(y)
    params = ((x.T*x).I*x.T*(y.T))
    b = params[0]
    w = params[1:]
    return w, b
    
    
if __name__ == "__main__":
    n_samples = 10000
    x_1, x_2 = np.random.random(n_samples),np.random.random(n_samples)
    theta_1, theta_2, b = map(float, input().split())
    
    x, y= np.vstack([x_1, x_2]), theta_1*x_1 + theta_2*x_2 + b
    params = getParams(x, y)

梯度下降逼近法

线性回归梯度下降法，这里的梯度下降是指计算损失函数(loss function)的梯度；这个损失函数通常是均方差函数
$$\begin{array}{lc}E(w)& =\frac{1}{2m}\sum _{i=1}^m({\hat{y}}^i-y^i{)}^2\\ & =\frac{1}{2m}\sum _{i=1}^m(w^T{x}^i+b-y^i{)}^2\end{array}$$
然后通过求解目标函数对于$w$的梯度和对于$b$的梯度，来不断的迭代更新他们的值，使得结果更加逼近真实值
$$\begin{gathered} \begin{array}{lc}\frac{\partial E}{\partial w}& =\frac{1}{m}\sum _{i=1}^m(w^T{x}^i+b-y^i)x^i\\ & =\frac{1}{m}\sum _{i=1}^m({\hat{y}}^i-y^i)x^i\end{array} \\ \begin{array}{lc}\frac{\partial E}{\partial b}& =\frac{1}{m}\sum _{i=1}^m(w^T{x}^i+b-y^i)x^i\\ & =\frac{1}{m}\sum _{i=1}^m({\hat{y}}^i-y^i)\end{array} \end{gathered}$$
迭代更新:
$$\begin{array}{l}w=w-\alpha \frac{\partial E}{\partial w}\\ \\ b=b-\alpha \frac{\partial E}{\partial b}\end{array}$$
其中$\alpha$表示学习率(learning rate)，用于控制每一步下降的步长(步长要选择适中，步长太小下降的速度比较慢，太大则可能造成不收敛；通常也可以选择动态步长的方法，即前期使用大步长用于快速下降，后期采用小步长使用收敛精度)，这里由于使用的是逼近的方法，所以我们可以设置通过循环的次数来结束求解过程，也可以使用设置阈值的方法结束求解。

代码实现

这里代码实现，分为使用numpy的方法和使用机器学习框架(pytorch)的方法

使用numpy

import numpy as np

def gradient_descent(X, y, lr=0.01, threshold=1e-3):
    y = y[:, np.newaxis]
    
    W =np.ones(shape=(1,X.shape[1]))
    b =np.array([[1]])
    n = y.shape[0]
    
    while True:
        pred_y = np.dot(X, W.T) + b
        loss = np.dot((y-pred_y).T, (y-pred_y))/(2*n)
        if np.abs(loss) < threshold:
            break
        
        w_gradient = -(1.0/n)*np.dot((y-pred_y).T, X)
        b_gradient = -(1.0/n)*np.dot((y-pred_y).T, np.ones(shape=[n, 1]))

        W = W - w_gradient
        b = b - b_gradient
        print(loss)
    
    return W[0][0], W[0][1], b[0][0]


    
if __name__ == "__main__":
    n_samples = 10000
    x_0 = np.ones(n_samples)
    x_1, x_2 = np.random.random(n_samples),np.random.random(n_samples)
    theta_1, theta_2, b = map(float, input().split())
    X, y= np.vstack([x_1, x_2]).T, theta_1*x_1 + theta_2*x_2 + b
    
    params = gradient_descent(X, y)
    result = []
    for param in params:
        result.append('{:.02f}'.format(param))
    print(' '.join(result))

使用pytorch

import torch
import torch.nn as nn
import numpy as np

# 线性回归
class LinerRegression(nn.Module):
    def __init__(self):
        super(LinerRegression, self).__init__()
        self.pred = torch.nn.Linear(2, 1)
    
    def forward(self, x):
        x = self.pred(x)
        return x

# 产生10000个数据
n_samples = 10000
x_1, x_2 = np.random.random(n_samples),np.random.random(n_samples)
x_train = np.vstack([x_1, x_2]).T
# 设置想要的权重，产生y
theta_1, theta_2, b = 8.69, -20.66, 76.12
y_train = theta_1*x_1 + theta_2*x_2 + b

# 把numpy array转为为tensor
x_train = torch.from_numpy(x_train)
y_train = torch.from_numpy(y_train)
# 为数据添加一个维度，和预测结果相匹配
y_train = y_train.unsqueeze(1)

# 调用网络
net = LinerRegression().double()

# 采用随机梯度下降法优化函数
optimizer = torch.optim.SGD(net.parameters(), lr = 0.01)
# 使用均方差作为 loss function
criterion = nn.MSELoss()

# 迭代循环，当loss小于一定的值时，结束
num = 0
while True:
    pred_y = net(x_train)
    loss = criterion(pred_y, y_train)
    
    optimizer.zero_grad()
    loss.backward()
    optimizer.step()
    
    if (num +1)%20==0:
        print('Loss: {:.6f}'.format(loss.data))
    
    if loss < 1e-3:
        break
    num += 1

# 保存参数
torch.save(net.state_dict(), 'params.pth')

# 加载参数
net.load_state_dict(torch.load('params.pth'))
# 输出每一个参数的值
for name, param in net.named_parameters():
    if name == 'pred.weight':
        print('W1:{:.02f}, W2:{:.02f}'.format(param.data[0][0].numpy(), param.data[0][1].numpy()))
    else:
        print('b:{:.02f}'.format(param.data[0].numpy()))

欢迎大家关注我的个人公众号，同样的也是和该博客账号一样，专注分享技术问题，我们一起学习进步

注： 文中有写错的地方，欢迎大家不吝指正！！!

作者：黑暗主宰

邮箱：shengzhanhe@gmail.com