利用最大似然估计和最小二乘法的直线拟合问题

1．最大似然估计原理

参数估计是给定观测的$n$个样本数据$(x_1,x_2,\cdots ,x_n)$，利用这些样本数据来估计未知的非随机常数参数$\theta$的问题。而最大似然估计则是参数估计问题中的一个非常常用的方法。其原理是所选择的$\theta$值使得观测样本数据发生的可能性最大。该可能性可以用联合概率密度$f(x_1,x_2,\cdots ,x_n;\theta )$（也称为似然函数）来表示，此时利用最大似然估计（MLE）所估计的参数可表示为
$$\left[{\hat{\theta }}_{ML}\right]=\underset{\theta }{\mathrm{argmax}}f(x_1,x_2,\cdots ,x_n;\theta )$$
一般对似然函数取对数会方便后续的求解，其对数似然函数为
$$L(x_1,x_2,\cdots ,x_n;\theta )\equiv \log f(x_1,x_2,\cdots ,x_n;\theta )$$
如果最大似然函数是可微并且存在上界，那么则可以通过求导的方式来求解参数${\hat{\theta }}_{ML}$，因为极值点的导数为零，即
$${\frac{\partial \log f(x_1,x_2,\cdots ,x_n;\theta )}{\partial \theta }|}_{\theta ={\hat{\theta }}_{ML}}=0$$

2．求解直线拟合问题

假设存在$n$个样本数据$(x_1,y_1),(x_2,y_2),\cdots ,(x_n,y_n)$，对其进行直线拟合使得
$$\begin{array}{cc}{y}_i=w{x}_i+b+{\epsilon }_i& i=1,2,\cdots \end{array}n$$
则利用最大似然估计所需估计的参数$\theta ={[w,b]}^T$，首先需要求取似然函数，假设残差${\epsilon }_i$服从均值为零的正态分布（高斯分布），那么${\epsilon }_i$的概率密度函数为
$$f({\epsilon }_i)=\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp (-\frac{{\epsilon }_i^2}{2{\sigma }^2})$$
由于 ${\epsilon }_i=y_i-w{x}_i-b$，因此
$$f(x_i,y_i;w,b)=\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp (-\frac{{(y_i-w{x}_i-b)}^2}{2{\sigma }^2})$$
假设各观测样本之间相互独立互不影响，那么$n$个样本数据之间的似然函数可以表示为
$$f(x_1,y_1,x_2,y_2,\cdots x_n,y_n,;w,b)=\prod _{i=1}^{n}f(x_i,y_i;w,b)$$
对数似然函数可表示为
$$\begin{array}{rl}& L(x_1,y_1,x_2,y_2,\cdots x_n,y_n,;w,b)=\log \prod _{i=1}^{n}f(x_i,y_i;w,b)\\ & =\log \prod _{i=1}^n\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp (-\frac{{(y_i-w{x}_i-b)}^2}{2{\sigma }^2})\\ & =\sum _{i=1}^n\log [\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp (-\frac{{(y_i-w{x}_i-b)}^2}{2{\sigma }^2})]\\ & =-\frac{n}{2}\log (2\pi {\sigma }^2)-\frac{1}{2{\sigma }^2}\sum _{i=1}^n{(y_i-w{x}_i-b)}^2\end{array}$$
可以看出标准差$\sigma$的取值不会影响参数$w,b$取何值时似然函数最大，因此使得上述对数似然函数最大可以简化为最小二乘问题，即
$$\min \sum _{i=1}^n{(y_i-w{x}_i-b)}^2$$
求解多元线性回归问题同样也具有相类似的推导，具体可参见https://zhuanlan.zhihu.com/p/143416436

3．最小二乘

定义向量$y={\left[y_1,y_2,\cdots ,y_n\right]}^T$，$\theta ={[w,b]}^T$，矩阵
$$A=\left[\begin{array}{cc}{x}_1& 1\\ x_2& 1\\ ⋮& ⋮\\ x_n& 1\end{array}\right]$$
则
$$\min \sum _{i=1}^n{(y_i-w{x}_i-b)}^2=\min {\Vert y-A\theta \Vert }^2$$
该最小二乘问题的解为求解如下法线方程组
$$A^{T}A\hat{\theta }=A^{T}y$$
$$\hat{\theta }={(A^{T}A)}^{-1}{A}^{T}y$$
所求得的参数$\hat{\theta }$即为直线拟合的最小二乘解，同时也是高斯分布假设下的最大似然估计，这两者是等效的。

4. 仿真

仿真参数设置为$y=3x+2+\epsilon$，$\epsilon$满足标准差$\sigma =0.5$的零均值正态分布，在$x\in (0,10)$的范围内随机选取$n=20$个点，所选取的点分布与直线拟合的结果如下

所求得的参数$\theta ={[3,0074,1.9817]}^T$。
代码如下：

clc;clear;close all
%% 设置仿真数据
n=20;%观测样本数
x=sort(rand(n,1))*10;
wucha=randn(n,1)*0.5;%生成满足方差为0.5的正态分布随机数
y=3*x+2+wucha;
figure;plot(x,y,'.','markersize',10);
ylim([0,35]);xlim([0,10]);
hold on;
%% 最小二乘法求解
A=[x,ones(n,1)];
thet=(A'*A)^-1*A'*y;
plot(x,A*thet,'-');

参考

《概率、随机变量与随机过程》第4版，帕普里斯等著；
https://zhuanlan.zhihu.com/p/143416436；
《线性代数》第9版，史蒂文J.利昂著