线性回归二三事

目录

一、线性回归定义

二、概率角度看线性回归

三、线性回归正则化

四、贝叶斯角度看正则化

五、线性回归扩展

六、从零实现线性回归代码分析

第一步，导入实验需要的包或模块

第二步，生成数据集

第三步、预览原始数据

第四步、数据读取

第五步、拟合参数

---

一、线性回归定义

 目标：给定一个数据集（具有样本和标签的有监督数据集），拟合出一条最符合给定数据的直线（本质就是拟合出一条最符合所有样本与其标签的函数，当然这里的直线是广义直线，数据样本一般都不是一维的）。

如上图，目标就是拟合出f(w)，即目标就是根据一个采样得到的数据集来估计出参数向量w的值。

参数估计方法：

最小二乘法，即参数估计中的最小二乘估计。

一种直观的解释是把每个样本点处的误差累加起来，最后使总误差最小。

二、概率角度看线性回归

最大似然估计：

点估计是频率派的做法，下面的最大似然估计就是最大似然点估计，从频率派的角度，参数w是一个未知但确定的值。所以式子里的分号代表以什么为参数，而不是条件概率。

似然是某个数据集的似然，顾名思义，似然，就是可能性，用最大似然估计的时候需要先已知（或假设出来）样本服从的概率分布，即需要已知概率密度函数，然后基于采样的iid假设，则我们采样得到整个样本集的可能性就是把每个样本代入概率密度函数，然后连乘。由于连乘之后求导求极值不方便，因此用对数函数将连乘变成连加。

 结论：

        当我们假设噪声分布（观测值与真实值的差的分布）是正态分布的时候，线性回归中参数向量的最小二乘估计就等价于最大似然估计。

三、线性回归正则化

正则化：

        就是在原损失函数上加一个关于参数的函数作为惩罚项，目的是解决过拟合。有时候甚至能使得解向量更稀疏。

        线性回归问题中一般有两种正则化方法，L1正则化和L2正则化。带有L1正则化的线性回归也被称为lasso，带有L2正则化的线性回归称为岭回归。

                                  

 lasso使得解向量稀疏，而岭回归使得模型一定有解析解。

 可以看到，如果样本数目小于维数的话，很有可能X^TX是不可逆的，但如果用了L2正则化，X^TX加上入I之后肯定是正定的，即可逆的，线性回归的解析解就一定存在了。

四、贝叶斯角度看正则化

         岭回归是最小化带有L2正则化的L2距离损失函数，它就等价于噪声服从高斯分布时候的最大后验。

噪声服从高斯分布：

        1：最小二乘估计等价于最大似然估计

        2：带有L2正则化的最小二乘估计等价于最大后验估计

        注意：所谓似然，这里其实是样本集的似然，也就是一个把样本集中的所有样本都代入这个样本集中样本服从的分布的密度函数，似然基于条件分布，这里的条件就是参数w（在贝叶斯角度w也是一个随机变量，所以可以作为条件，频率派角度就不行）。因为假设样本是独立同分布的，所以可以把所有样本代入密度函数再连乘，然后加入log函数不改变关于w的函数的单调性，把连乘变成连加。最大后验估计就是在样本x服从的条件分布（条件就是参数w）的密度函数的基础上，再乘以一个条件也就是参数w服从的先验分布。构成后验表达式。因为我们做最小二成估计也好，最大似然估计也好，最大后验估计也好，最终目的都是估计参数w的值，所以我们提前假设一个参数w的分布，当然可以看作是w的先验，因为在我们解出来w的值之前，就已经知道关于w的信息了，所以我们用L2正则化的时候，对w进行约束，也相当于是w的一个先验。

        L1正则化相比于L2更容易获得稀疏解，使得w有更少的非零分量，从而间接完成了降维，因为w中的0对应的那些特征都不会用上，相当于对特征进行降维了。

五、线性回归扩展

对数线性回归：

        如果样本的输出是在指数尺度上变化的，就可以用对数线性回归。

        lny = w^Tx + b，即 y = exp{w^Tx + b}

        为什么叫对数线性回归？因为输出的对数lny与输入x之间是线性关系，但y本身是在x的指数尺度上变化的。这样的模型也叫广义线性模型。

        上边exp{}指数函数也叫激活函数，即对原先关于x的线性输出做一次非线性变换，而从线性回归的角度来说，ln函数被称为link function，虽然它们两个是同一个东西，只不过放在等号的左边或者右边，但是从不同的角度有不同的叫法。

对数几率回归：

        在对数线性回归中，如果用sigmoid函数（也叫对数几率函数）作为激活函数，此时对数线性回归称作对数几率回归。其实所有s形函数都叫sigmoid函数对数几率函数只是最具代表性的一种）

        y = 1/{1+exp{-(w^Tx +b)}}

        即把w^Tx +b代入到sigmoid函数的变量z中。

        此时sigmoid函数把输出值变换到【0，1】区间，并且输出值在0.5附近变换很陡。

        把对数几率回归变换成link函数表示的形式，我们就能看到其名字的由来：

        ln(y / (1 - y)) = w^Tx + b

        因为y的值域是【0，1】，可以把y看成样本x是正例的可能性，1 - y 则是样本是反例的可能性。正例可能性与反例可能性的笔试 y / (1 - y) 被称为几率 odds，反应了样本x是正例的相对可能性，再取对数ln，则称为对数几率。因为对数几率回归要回归的东西就是对数几率。

特点：

1：直接对分类可能性进行建模，无需实现假设数据分布，避免了假设分布不准确带来的问题。 2：不仅能预测类别，而是得到近似概率预测

3：对率函数是任意阶可导的凸函数，数学性质好。

六、从零实现线性回归代码分析

https://tangshusen.me/Dive-into-DL-PyTorch/#/chapter03_DL-basics/3.2_linear-regression-scratch

第一步，导入实验需要的包或模块

%matplotlib inline
import torch
from IPython import display
from matplotlib import pyplot as plt
import numpy as np
import random

原实验是在Jupyter notebook中做的，其本质是一个python的交互式shell

%matplotlib inline是IPython的魔法命令（https://blog.youkuaiyun.com/liming89/article/details/109662966），其作用是当用matplotlib画图的时候可以直接内嵌显示，不用生成额外的图像，这样在用plt画图的时候就不用plt.show()来手动画出来了，IPython内嵌了许多的魔法命令。

torch即著名机器学习库。

display是IPython中用来设置显示参数的API，本质是一个.py模块。是IPython独有的。从它的位置也能看出来，site-packages目录下一般用于放python的库，而display.py是在IPython目录下的。

matplotlib是一个python绘图库，通常和numpy一起使用可以作为matlib的替代方案，pyplot是matplotlib中的绘图.py模块。

numpy是python的一个数学函数库，可用于各种矩阵运算。

random是一个.py标准库模块，用于生成随机数以及将列表中的元素打乱等等。

第二步，生成数据集

# 生成数据集
num_inputs = 2 # 特征向量维度
num_examples = 1000 # 训练样本数
true_w = [2, -3.4] # 真实系数
true_b = 4.2 # 真实偏置项
features = torch.randn(num_examples, num_inputs, dtype=torch.float32)
labels = true_w[0] * features[:, 0] + true_w[1] * features[:, 1] + true_b
labels += torch.from_numpy(np.random.normal(0, 0.01, size=labels.size()))
print(labels)

torch.randn()函数用于生成某种维度的服从正态分布的张量。这里生成了1000行2列的数据，每一行是一个样本特征向量，其特征维度是2。

labels也是一个张量，注意python在运算的时候有广播机制，因此张量也可以直接与标量相加。

在生成标签的时候，先用numpy的random模块中的normal函数，生成了均值为0，标准差为0.01，size与labels相同的数组，然后转换成torch中的tensor，相当于给标签加了噪声。这里需要注意，python标准库中有一个random模块，一般不用于生成序列类随机数，random模块中的方法多用于生成一个随机数，或者打乱列表类的序列。torch中也有一个torch.random，但是也不用于生成某种分布的张量，用于生成张量的一般是torch.randn这类的函数。但是torch.randn是生成标准分布的，所以要生成确定形状，从确定方差和均值的正太分布中采样得到的张量还是需要先用numpy.random.normal来生成一个多维数组，然后再转换成张量。

第三步、预览原始数据

def use_svg_display():
    display.set_matplotlib_formats('svg') # 该函数应该就是设置格式，设置成用矢量图表示
def set_figsize(figsize=(3.5, 2.5)):
    use_svg_display() # 顾名思义，用矢量图显示
    plt.rcParams['figure.figsize'] = figsize # 设置图的尺寸
set_figsize()
#fig = plt.gcf()
plt.scatter(features[:, 0].numpy(), labels.numpy()) # sactter的两个参数分别是x和y，即散点图的横纵坐标
plt.scatter(features[:, 1].numpy(), labels.numpy())
plt.savefig(r"D:\常用\notebook\fig.svg")

display是IPython中用于显示图像设置的模块，这里display.set_matplotlib_formats是设置图像的格式，这里设置成用矢量图显示。

plt.rcParams是一个字典对象，里边有很多关于画布的设置，其中画布尺寸就是这里的figure.figsize，这个不用管也行，默认就是（3.5，2.5）

从上图也可以看出，第一个标签关于第一个特征大概是正相关，关于第二个特征是负相关，这与我们生成标签时候用的系数（2，-3.4）也是符合的。

第四步、数据读取

# 数据读取函数
def data_iter(batch_size, features, labels):
    num_examples = len(features) # 样本数
    indices = list(range(num_examples)) # range返回一个可迭代对象，然后用list将该可迭代对象转化成list
    random.shuffle(indices)  # 参数可以是列表，作用是将一个序列的所有元素打乱顺序
    for i in range(0, num_examples, batch_size):
        j = torch.LongTensor(indices[i : min(i + batch_size, num_examples)]) # i 是一个batch的起始序号，min()是结束序号，最后一个可能不足一个batch
        # 因为indices是打乱过的，因此j是一个长度为batch_size的1阶张量，可以认为是features张量的索引，用来从features中抽取随机的batch_size个数据
        # 因为indices是从0到num_examples的整数列表，里边没有重复的数，因此我们的batch之间也是不重叠的。
        yield features.index_select(0, j), labels.index_select(0, j) # index_select中第一个参数0表示按行索引，j是一个张量，表示用来索引的索引号
        # yield把函数改造成了一个生成器，相当于data_iter函数在for循环语句中买个返回一个张量，即索引得到的张量

 num_examples是样本的数量，range（）是python自带的函数，用于生成一个可迭代对象，默认是[0, num_examples)即左闭右开区间的一个可迭代对象。然后我们用list()把这个可迭代对象转换成一个列表indices，其中每个元素都是一个索引相当于。

再用python标准库中的random模块将列表indices打乱。

接下来for语句，range生成可迭代对象，从0开始，步长是batch_size，最后是num_examples - 1。i和min()用来从indices中每次切片batch_size个元素，因此得到的j是大小为batch_size的一个张量（最后一次循环的j可能小于batch_size）

然后用index_select函数来从张量中按索引来选取出子张量，其中第一个参数0代表按行索引，第二个参数j就是索引用的列表。

注意，看到yield就可以知道，函数data_iter已经被改造成了生成器，生成器是迭代器的一种特殊情况，迭代器是可迭代对象的一种特殊情况。

batch_size = 10 # 读取第一个batch的数据看看
for X, y in data_iter(batch_size, features, labels):
    print(X, y)
    break

因为data_iter是一个生成器，因此我们就可以用for in 的方式来遍历了。

第五步、拟合参数

w = torch.tensor(np.random.normal(0, 0.01, (num_inputs, 1)), dtype=torch.float32)
b = torch.zeros(1, dtype=torch.float32)
w.requires_grad_(requires_grad=True)
b.requires_grad_(requires_grad=True)

w是系数张量，形状是（2，1），因为样本有两个特征

b是偏置项。

张量w和b都是需要用梯度下降来更新的，因此其成员requires_grad都设置为true。

def linreg(X, w, b):
    return torch.mm(X, w) + b
def squared_loss(y_hat, y):
    return (y_hat - y.view(y_hat.size())) ** 2 / 2

X形状（batch_size, 2）, w形状（2，1），torch.mm是矩阵乘法（非点乘）。

因此linreg函数就是相当于前向传播，返回的结果即预测值。

squared_loss即平方差损失函数。

def sgd(params, lr, batch_size):
    for param in params:
        print(param)
        param.data -= lr * param.grad / batch_size # params就是w和b，注意这里用了param.data，也就是只改变大小，不计入计算图

这是参数优化方法sgd，params中有两个张量，一个是张量w，一个是张量b，其中param.data与param.grad是形状相同的张量。只所以用param.data而不是直接用张量操作，是为了在优化参数的时候防止影响该参数张量对应的梯度。

lr = 0.03
num_epochs = 3
net = linreg # python中可以直接把函数作为一个对象来看待
loss = squared_loss

for epoch in range(num_epochs):
    for X, y in data_iter(batch_size, features, labels):
        l = loss(net(X, w, b), y).sum()
        l.backward()
        sgd([w, b], lr, batch_size) # 上一步计算出梯度，这一步用minisgd来进行参数更新
        w.grad.data.zero_()
        b.grad.data.zero_() # 更新完之后将本轮梯度清零
    train_l = loss(net(features, w, b), labels)
    print('epoch %d, loss %f' % (epoch + 1, train_l.mean().item()))

需要注意每次循环之后梯度清零。