利用最大似然估计和最小二乘法的直线拟合问题
即为直线拟合的最小二乘解,同时也是高斯分布假设下的最大似然估计,这两者是等效的。而最大似然估计则是参数估计问题中的一个非常常用的方法。(也称为似然函数)来表示,此时利用最大似然估计(MLE)所估计的参数可表示为。取何值时似然函数最大,因此使得上述对数似然函数最大可以简化为最小二乘问题,即。如果最大似然函数是可微并且存在上界,那么则可以通过求导的方式来求解参数。一般对似然函数取对数会方便后续的求解,
利用最大似然估计和最小二乘法的直线拟合问题
1.最大似然估计原理
参数估计是给定观测的 n n n个样本数据 ( x 1 , x 2 , ⋯ , x n ) ({{x}_{1}},{{x}_{2}},\cdots ,{{x}_{n}}) (x1,x2,⋯,xn),利用这些样本数据来估计未知的非随机常数参数 θ \theta θ的问题。而最大似然估计则是参数估计问题中的一个非常常用的方法。其原理是所选择的 θ \theta θ值使得观测样本数据发生的可能性最大。该可能性可以用联合概率密度 f ( x 1 , x 2 , ⋯ , x n ; θ ) f({{x}_{1}},{{x}_{2}},\cdots ,{{x}_{n}};\theta ) f(x1,x2,⋯,xn;θ)(也称为似然函数)来表示,此时利用最大似然估计(MLE)所估计的参数可表示为
[ θ ^ M L ] = arg max θ f ( x 1 , x 2 , ⋯ , x n ; θ ) \left[ {{{\hat{\theta }}}_{ML}} \right]=\underset{\theta }{\mathop{\arg \max }}\,f({{x}_{1}},{{x}_{2}},\cdots ,{{x}_{n}};\theta )\ [θ^ML]=θargmaxf(x1,x2,⋯,xn;θ)
一般对似然函数取对数会方便后续的求解,其对数似然函数为
L ( x 1 , x 2 , ⋯ , x n ; θ ) ≡ log f ( x 1 , x 2 , ⋯ , x n ; θ ) L({{x}_{1}},{{x}_{2}},\cdots ,{{x}_{n}};\theta )\equiv \log f({{x}_{1}},{{x}_{2}},\cdots ,{{x}_{n}};\theta ) L(x1,x2,⋯,xn;θ)≡logf(x1,x2,⋯,xn;θ)
如果最大似然函数是可微并且存在上界,那么则可以通过求导的方式来求解参数 θ ^ M L {{\hat{\theta }}_{ML}} θ^ML,因为极值点的导数为零,即
∂ log f ( x 1 , x 2 , ⋯ , x n ; θ ) ∂ θ ∣ θ = θ ^ M L = 0 {{\left. \frac{\partial \log f({{x}_{1}},{{x}_{2}},\cdots ,{{x}_{n}};\theta )}{\partial \theta } \right|}_{\theta ={{{\hat{\theta }}}_{ML}}}}=0 ∂θ∂logf(x1,x2,⋯,xn;θ)
θ=θ^ML=0
2.求解直线拟合问题
假设存在 n n n个样本数据 ( x 1 , y 1 ) , ( x 2 , y 2 ) , ⋯ , ( x n , y n ) ({{x}_{1}},{{y}_{1}}),({{x}_{2}},{{y}_{2}}),\cdots ,({{x}_{n}},{{y}_{n}}) (x1,y1),(x2,y2),⋯,(xn,yn),对其进行直线拟合使得
y i = w x i + b + ε i i = 1 , 2 , ⋯ n \begin{matrix} {{y}_{i}}=w{{x}_{i}}+b+{{\varepsilon }_{i}} & i=1,2,\cdots \\ \end{matrix}n yi=wxi+b+εii=1,2,⋯n
则利用最大似然估计所需估计的参数 θ = [ w , b ] T \theta ={{[w,b]}^{T}} θ=[w,b]T,首先需要求取似然函数,假设残差 ε i {{\varepsilon }_{i}} εi服从均值为零的正态分布(高斯分布),那么 ε i {{\varepsilon }_{i}} εi的概率密度函数为
f ( ε i ) = 1 2 π σ 2 exp ( − ε i 2 2 σ 2 ) f({{\varepsilon }_{i}})=\frac{1}{\sqrt{2\pi {{\sigma }^{2}}}}\exp (-\frac{\varepsilon _{i}^{2}}{2{{\sigma }^{2}}}) f(εi)=2πσ21exp(−2σ2εi2)
由于 ε i = y i − w x i − b {{\varepsilon }_{i}}={{y}_{i}}-w{{x}_{i}}-b εi=yi−wxi−b,因此
f ( x i , y i ; w , b ) = 1 2 π σ 2 exp ( − ( y i − w x i − b ) 2 2 σ 2 ) 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}}}) f(xi,yi;w,b)=2πσ21exp(−2σ2(yi−wxi−b)2)
假设各观测样本之间相互独立互不影响,那么 n n n个样本数据之间的似然函数可以表示为
f ( x 1 , y 1 , x 2 , y 2 , ⋯ x n , y n , ; w , b ) = ∏ i = 1 n f ( x i , y i ; w , b ) f({{x}_{1}},{{y}_{1}},{{x}_{2}},{{y}_{2}},\cdots {{x}_{n}},{{y}_{n}},;w,b)=\prod\limits_{i=1}^{n}{f({{x}_{i}},{{y}_{i}};w,b)} f(x1,y1,x2,y2,⋯xn,yn,;w,b)=i=1∏nf(xi,yi;w,b)
对数似然函数可表示为
L ( x 1 , y 1 , x 2 , y 2 , ⋯ x n , y n , ; w , b ) = log ∏ i = 1 n f ( x i , y i ; w , b ) = log ∏ i = 1 n 1 2 π σ 2 exp ( − ( y i − w x i − b ) 2 2 σ 2 ) = ∑ i = 1 n log [ 1 2 π σ 2 exp ( − ( y i − w x i − b ) 2 2 σ 2 ) ] = − n 2 log ( 2 π σ 2 ) − 1 2 σ 2 ∑ i = 1 n ( y i − w x i − b ) 2 \begin{align} & L({{x}_{1}},{{y}_{1}},{{x}_{2}},{{y}_{2}},\cdots {{x}_{n}},{{y}_{n}},;w,b)=\log \prod\limits_{i=1}^{n}{f({{x}_{i}},{{y}_{i}};w,b)} \\ & =\log \prod\limits_{i=1}^{n}{\frac{1}{\sqrt{2\pi {{\sigma }^{2}}}}\exp (-\frac{{{({{y}_{i}}-w{{x}_{i}}-b)}^{2}}}{2{{\sigma }^{2}}})} \\ & =\sum\limits_{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\limits_{i=1}^{n}{{{({{y}_{i}}-w{{x}_{i}}-b)}^{2}}} \\ \end{align} L(x1,y1,x2,y2,⋯xn,yn,;w,b)=logi=1∏nf(xi,yi;w,b)=logi=1∏n2πσ21exp(−2σ2(yi−wxi−b)2)=i=1∑nlog[2πσ21exp(−2σ2(yi−wxi−b)2)]=−2nlog(2πσ2)−2σ21i=1∑n(yi−wxi−b)2
可以看出标准差 σ \sigma σ的取值不会影响参数 w , b w,b w,b取何值时似然函数最大,因此使得上述对数似然函数最大可以简化为最小二乘问题,即
min ∑ i = 1 n ( y i − w x i − b ) 2 \min \sum\limits_{i=1}^{n}{{{({{y}_{i}}-w{{x}_{i}}-b)}^{2}}} mini=1∑n(yi−wxi−b)2
求解多元线性回归问题同样也具有相类似的推导,具体可参见https://zhuanlan.zhihu.com/p/143416436
3.最小二乘
定义向量 y = [ y 1 , y 2 , ⋯ , y n ] T y={{\left[ {{y}_{1}},{{y}_{2}},\cdots ,{{y}_{n}} \right]}^{T}} y=[y1,y2,⋯,yn]T, θ = [ w , b ] T \theta ={{[w,b]}^{T}} θ=[w,b]T,矩阵
A = [ x 1 1 x 2 1 ⋮ ⋮ x n 1 ] A=\left[ \begin{matrix} {{x}_{1}} & 1 \\ {{x}_{2}} & 1 \\ \vdots & \vdots \\ {{x}_{n}} & 1 \\ \end{matrix} \right] A=
x1x2⋮xn11⋮1
则
min ∑ i = 1 n ( y i − w x i − b ) 2 = min ∥ y − A θ ∥ 2 \min \sum\limits_{i=1}^{n}{{{({{y}_{i}}-w{{x}_{i}}-b)}^{2}}}=\min {{\left\| y-A\theta \right\|}^{2}} mini=1∑n(yi−wxi−b)2=min∥y−Aθ∥2
该最小二乘问题的解为求解如下法线方程组
A T A θ ^ = A T y {{A}^{T}}A\hat{\theta }={{A}^{T}}y ATAθ^=ATy
θ ^ = ( A T A ) − 1 A T y \hat{\theta }={{({{A}^{T}}A)}^{-1}}{{A}^{T}}y θ^=(ATA)−1ATy
所求得的参数 θ ^ \hat{\theta } θ^即为直线拟合的最小二乘解,同时也是高斯分布假设下的最大似然估计,这两者是等效的。
4. 仿真
仿真参数设置为 y = 3 x + 2 + ε y=3x+2+\varepsilon y=3x+2+ε, ε \varepsilon ε满足标准差 σ = 0.5 \sigma =0.5 σ=0.5的零均值正态分布,在 x ∈ ( 0 , 10 ) x\in (0,10) x∈(0,10)的范围内随机选取 n = 20 n=20 n=20个点,所选取的点分布与直线拟合的结果如下
所求得的参数 θ = [ 3 , 0074 , 1.9817 ] T \theta ={{[3,0074,1.9817]}^{T}} θ=[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.利昂著
更多推荐
所有评论(0)