【数理知识】最小二乘法，从线性回归出发，数值举例并用最小二乘法求解回归模型

序号	内容
1	【数理知识】自由度 degree of freedom 及自由度的计算方法
2	【数理知识】刚体 rigid body 及刚体的运动
3	【数理知识】刚体基本运动，平动，转动
4	【数理知识】向量数乘，内积，外积，matlab代码实现
5	【数理知识】协方差，随机变量的的协方差，随机变量分别是单个数字和向量时的协方差
6	【数理知识】旋转矩阵的推导过程，基于向量的旋转来实现，同时解决欧式变换的非线性局限

查了下自己的笔记，之前在用 SG 滤波器时提到过最小二乘法，这里放一下文章链接方便自己查阅：【UWB】Savitzky Golay filter SG滤波器原理讲解。

@[toc]

最小二乘法（Least Squared Method）是一种数学优化技术，用于最小化预测值与真实值之间的差异（通常表示为残差平方和的形式）。

最小二乘法的核心思想是：通过最小化预测值和真实值之间的差异的平方和，来找到模型最佳的参数值。

因为最小二乘法就是一门优化技术，单纯的讨论最小二乘法比较枯燥。因此，将最小二乘法应用在回归分析中的线性回归，帮助求解出回归模型的最佳参数，这样的方式可以更好的理解最小二乘法的思想。

回归分析的最初目的是估计模型的参数以便达到对数据的最佳拟合。在很多情况下，特别是在线性回归中，最小二乘法可以为我们提供一个封闭的、解析的解决方案来寻找最佳参数。因此，最小二乘法经常用于拟合直线和多项式得到一组数据点，从而预测或解释变量的行为。

接下来，将先阐述线性回归，并阐述两个简单的线性回归模型。然后举具体的例子，通过最小二乘法的思想求解出最佳模型参数，从而实现对使用最小二乘法求解模型参数的理解。

1. 线性回归

讨论线性回归（Linear Regression）时，同时结合回归分析（Regression Analysis）一起讨论下。

线性回归和回归分析都是统计学中常用的方法，但它们之间存在一些关键的区别：

在定义上，

线性回归：线性回归是一个预测方法，它试图找到一个线性函数（对于简单线性回归）或多个线性函数（对于多元线性回归），这些函数可以最好地描述两个（或更多）变量之间的关系。
回归分析：回归分析是一个更广泛的术语，用于描述一个变量（或多个变量）与一个或多个其他变量之间关系的统计方法。线性回归只是回归分析中的一种形式。

在类型方面，

线性回归：主要关注线性关系。
回归分析：可以包括线性回归、多项式回归、逻辑回归、岭回归等多种类型。

在目的上，

线性回归：预测或解释一个响应变量和一个或多个预测变量之间的关系。
回归分析：探索和建模变量之间的关系，这可能是线性的、非线性的，或者是其他的关系。

在应用上，

线性回归：当我们认为变量之间存在线性关系时，通常使用线性回归。
回归分析：可以应用于各种关系，包括但不限于线性关系。

简而言之，线性回归是回归分析的一个子集。回归分析包括多种方法，用于模拟和解释变量之间的关系，而线性回归专门关注线性关系。

给一个随机样本（ $y, x_1, x_2, \cdots$ ），一个线性回归模型假设
响应变量（通常称为因变量或目标） $y$ 和
解释变量（通常称为自变量或特征） $x_1, x_2, \cdots$
之间的关系除了受到解释变量的影响之外，还有其它变量的存在。我们加入一个误差项 $\epsilon$ （也是一个随机变量）来捕获除了 $x_1, x_2, \cdots$ 之外任何对 $y$ 的影响。

不过在本文的分析中，我们主要是为了学习回归模型，故简化掉误差项 $\epsilon$ 的影响。

同时，根据解释变量数目的不同，线性回归模型还分为了简单线性回归和多元线性回归。

1. 简单线性回归模型

描述了一个响应变量 $y$ 和一个解释变量 $x$ 之间的关系。

假设这种关系是线性的，也就是一次函数

$\begin{aligned} y &= \beta _0 + \beta_1 x \end{aligned}$

其中 $y$ 是响应变量， $x$ 是解释变量， $\beta_0$ 和 $\beta_1$ 是回归系数。

2. 多元线性回归模型

描述了一个响应变量 $y$ 和两个或两个以上解释变量 $x_1, x_2, \cdots$ 之间的关系。

$\begin{aligned} y &= \beta _0 + \beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_p x_p \end{aligned}$

2. 从线性回归到最小二乘法，数值举例并用最小二乘法求解回归模型

最小二乘法的核心思想是寻找能使预测误差平方和最小化的模型参数。

具体来说，考虑一个简单的线性模型：

$\begin{aligned} y &= \beta _0 + \beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_p x_p \end{aligned}$

我们的目标是找到一组 $\beta$ 值，使得预测值 $\hat{y} = \beta _0 + \beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_p x_p$ 与实际值 $y$ 之间的差异最小。这种差异可以通过残差平方和（RSS）来衡量：

$\text{RSS} = \sum_{i=1}^{n} (y - \hat{y}_i)^2$

最小二乘法的目标是最小化 RSS。具体来说就是通过寻找能最小化 RSS 的 $\beta$ 值来实现。

最小二乘法背后的直观思想是尽量找到一个模型，使得预测值尽量接近观测值。这种方法是在基于误差平方损失的情况下得出的。

接下来通过举具体的数值例子，观察最小二乘法的作用。

1. 单解释变量，单响应变量

这种单解释变量单响应变量的回归模型，可以用公式表征为

$y = \beta_0 + \beta_1 x$

假设一个简单线性回归模型 $y = \beta_0 + \beta_1 x$ ，在实验中得到了四个数据点（ $x,y$ ）：（ $1,6$ ）（ $2,5$ ）（ $3,7$ ）（ $4,10$ ），找到最佳的匹配参数 $\beta_0, \beta_1$ 。

分别将数据点代入到此线性回归模型中有
$\begin{aligned} \beta_0 + \beta_1 1 &= 6 \\ \beta_0 + \beta_1 2 &= 5 \\ \beta_0 + \beta_1 3 &= 7 \\ \beta_0 + \beta_1 4 &= 10 \end{aligned}$

最小二乘法采用的方法就是尽量使等号两边差的平方最小，也就是找出如下函数的最小值：
$\begin{aligned} S(\beta_0, \beta_1) &= [6 - (\beta_0 + \beta_1 1)]^2 \\ &+ [5 - (\beta_0 + \beta_1 2)]^2 \\ &+ [7 - (\beta_0 + \beta_1 3)]^2 \\ &+ [10 - (\beta_0 + \beta_1 4)]^2 \end{aligned}$

% 初始化符号变量
syms beta_0 beta_1 real

% 定义你的公式
f = (6 - (beta_0 + beta_1 * 1))^2 ... 
  + (5 - (beta_0 + beta_1 * 2))^2 ... 
  + (7 - (beta_0 + beta_1 * 3))^2 ... 
  + (10 - (beta_0 + beta_1 * 4))^2;

% 展开公式
f_expanded = expand(f);

% 你可以尝试进一步简化它，但是否能简化到所需的形式是不确定的
f_simplified = simplify(f_expanded);

% 输出结果
disp(f_expanded);
disp(f_simplified);

>>
4*beta_0^2 + 20*beta_0*beta_1 - 56*beta_0 + 30*beta_1^2 - 154*beta_1 + 210
 
4*beta_0^2 + 20*beta_0*beta_1 - 56*beta_0 + 30*beta_1^2 - 154*beta_1 + 210

最小值可以通过对 $S(\beta_0, \beta_1)$ 分别求 $\beta_0$ 和 $\beta_1$ 的偏导数，并令其等于零得到。

$\begin{aligned} \frac{\partial S(\beta_0, \beta_1)}{\partial \beta_0} &= 8 \beta_0 + 20 \beta_1 - 56 = 0 \\ \frac{\partial S(\beta_0, \beta_1)}{\partial \beta_1} &= 20 \beta_0 + 60 \beta_1 - 154 = 0 \end{aligned}$

% 对 b0 求偏导数
df_dbeta_0 = diff(f, beta_0);

% 对 b1 求偏导数
df_dbeta_1 = diff(f, beta_1);

disp(df_dbeta_0)

disp(df_dbeta_1)

>> 
8*beta_0 + 20*beta_1 - 56
20*beta_0 + 60*beta_1 - 154

通过求解上述二元一次方程可以得到

$\begin{aligned} \beta_0 &= 3.5 \\ \beta_1 &= 1.4 \end{aligned}$

% 设置偏导数等于零并求解
solutions = solve([df_dbeta_0 == 0, df_dbeta_1 == 0], [beta_0, beta_1]);

% 输出解
disp(solutions.beta_0);

disp(solutions.beta_1);

>>
7/2
 
7/5

因此最佳回归模型为 $y = 3.5 + 1.4 x$ 。

x = [1  2  3  4]';
y = [6  5  7  10]';

plot_x = linspace(1,4,50);
plot_y = 3.5 + 1.4 * plot_x;

scatter(x(1), y(1)); hold on;
scatter(x(2), y(2));
scatter(x(3), y(3));
scatter(x(4), y(4));
plot(plot_x, plot_y);
xlabel("$x$", "Interpreter","latex", "FontSize",16);
ylabel("$y$", "Interpreter","latex", "FontSize",16);
grid on;

image.png

2. 多解释变量，单响应变量

这种多解释变量单响应变量的回归模型，可以用公式表征为

$y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_p x_p$

同样基于上述回归模型，我们这里取解释变量的个数为 $p = 2$ ，那么有 $y = \beta_0 + \beta_1 x_1 + \beta_2 x_2$ ，在实验中得到了四个数据点（ $x_1, x_2, y$ ）：（ $1,3,6$ ）（ $2,4,5$ ）（ $3,5,7$ ）（ $4,8,10$ ），找到最佳的匹配参数 $\beta_0, \beta_1, \beta_2$ 。

分别将数据点代入到此线性回归模型中有
$\begin{aligned} \beta_0 + \beta_1 1 + \beta_2 3 &= 6 \\ \beta_0 + \beta_1 2 + \beta_2 4&= 5 \\ \beta_0 + \beta_1 3 + \beta_2 5&= 7 \\ \beta_0 + \beta_1 4 + \beta_2 8&= 10 \end{aligned}$

得到差的函数有

$\begin{aligned} S(\beta_0, \beta_1, \beta_2) &= [6 - (\beta_0 + \beta_1 1 + \beta_2 3)]^2 \\ &+ [5 - (\beta_0 + \beta_1 2 + \beta_2 4)]^2 \\ &+ [7 - (\beta_0 + \beta_1 3 + \beta_2 5)]^2 \\ &+ [10 - (\beta_0 + \beta_1 4 + \beta_2 8)]^2 \end{aligned}$

syms beta_0 beta_1 beta_2 real

f = (6 - (beta_0 + beta_1 * 1 + beta_2 * 3))^2 ...
  + (5 - (beta_0 + beta_1 * 2 + beta_2 * 4))^2 ...
  + (7 - (beta_0 + beta_1 * 3 + beta_2 * 5))^2 ... 
  + (10 - (beta_0 + beta_1 * 4 + beta_2 * 8))^2;

f_expanded = expand(f);

disp(f_expanded);

>> 
4*beta_0^2 + 20*beta_0*beta_1 + 40*beta_0*beta_2 - 56*beta_0 + 30*beta_1^2 + 116*beta_1*beta_2 - 154*beta_1 + 114*beta_2^2 - 306*beta_2 + 210

分别求偏导数，并令其等于零
$\begin{aligned} \frac{\partial S(\beta_0, \beta_1, \beta_2)}{\partial \beta_0} &= 8 \beta_0 + 20 \beta_1 + 40 \beta_2 - 56 = 0 \\ \frac{\partial S(\beta_0, \beta_1, \beta_2)}{\partial \beta_1} &= 20 \beta_0 + 60 \beta_1 + 116 \beta_2 - 154 = 0 \\ \frac{\partial S(\beta_0, \beta_1, \beta_2)}{\partial \beta_2} &= 40 \beta_0 + 116 \beta_1 + 228 \beta_2 - 306 = 0 \end{aligned}$

df_dbeta_0 = diff(f, beta_0);
df_dbeta_1 = diff(f, beta_1);
df_dbeta_2 = diff(f, beta_2);

disp(df_dbeta_0)
disp(df_dbeta_1)
disp(df_dbeta_2)

>> 
8*beta_0 + 20*beta_1 + 40*beta_2 - 56
20*beta_0 + 60*beta_1 + 116*beta_2 - 154
40*beta_0 + 116*beta_1 + 228*beta_2 - 306

求解上述三元一次方程可以得到

$\begin{aligned} \beta_0 &= 2 \\ \beta_1 &= -1 \\ \beta_2 &= 1.5 \end{aligned}$

solutions = solve([df_dbeta_0 == 0, df_dbeta_1 == 0, df_dbeta_2 == 0], [beta_0, beta_1, beta_2]);

disp(solutions.beta_0);
disp(solutions.beta_1);
disp(solutions.beta_2);

>>
2
-1
3/2

因此最佳回归模型为 $y = 2 - x_1 + 1.5 x_2$ 。

x_1 = [1  2  3  4]';
x_2 = [3  4  5  8]';
y = [6  5  7  10]';

plot_x_1 = linspace(1,4,50);
plot_x_2 = linspace(3,8,50);
plot_y = 2 - 1 * plot_x_1 + 1.5 * plot_x_2;

figure()
subplot(2,2,1)
scatter3(x_1(1), x_2(1), y(1)); hold on;
scatter3(x_1(2), x_2(2), y(2));
scatter3(x_1(3), x_2(3), y(3));
scatter3(x_1(4), x_2(4), y(4));
plot3(plot_x_1, plot_x_2, plot_y);
xlabel("$x_1$", "Interpreter","latex", "FontSize",16);
ylabel("$x_2$", "Interpreter","latex", "FontSize",16);
zlabel("$y$", "Interpreter","latex", "FontSize",16);
grid on;

subplot(2,2,2)
scatter(x_1(1), y(1)); hold on;
scatter(x_1(2), y(2));
scatter(x_1(3), y(3));
scatter(x_1(4), y(4));
plot(plot_x_1, plot_y);
xlabel("$x_1$", "Interpreter","latex", "FontSize",16);
ylabel("$y$", "Interpreter","latex", "FontSize",16);
grid on;

subplot(2,2,3)
scatter(x_2(1), y(1)); hold on;
scatter(x_2(2), y(2));
scatter(x_2(3), y(3));
scatter(x_2(4), y(4));
plot(plot_x_2, plot_y);
xlabel("$x_2$", "Interpreter","latex", "FontSize",16);
ylabel("$y$", "Interpreter","latex", "FontSize",16);
grid on;

subplot(2,2,4)
scatter(x_1(1), x_2(1)); hold on;
scatter(x_1(2), x_2(2));
scatter(x_1(3), x_2(3));
scatter(x_1(4), x_2(4));
plot(plot_x_1, plot_x_2);
xlabel("$x_1$", "Interpreter","latex", "FontSize",16);
ylabel("$x_2$", "Interpreter","latex", "FontSize",16);
grid on;