Matlab实现最小二乘法的高斯牛顿迭代算法

最新推荐文章于 2023-12-19 21:51:33 发布

Antash

最新推荐文章于 2023-12-19 21:51:33 发布

阅读量1.6k

点赞数 3

文章标签： matlab 算法最小二乘法

本文链接：https://blog.csdn.net/Andrew1510/article/details/129961609

版权

一、算法解释

给定一组观测值（测量值） $(x_{i},y_{i})$ ，以及观测方程 $Y=H(x,\theta)+v$ ，求 $\theta$ .

对任意给定的初值 $\theta_{0}$ ，观测方程在 $\theta_0$ 处展开有

$Y=H(\theta)+\frac{\partial H}{\partial\theta}|_{\theta=\theta_{0}}(\theta - \theta_{0})+o^{1}+v$

也就是

$Y\approx H(\theta)+J(\theta_{0})(\theta - \theta_{0})+v$

忽略高阶小量并记 $\theta-\theta_{0}=\Delta \theta; Y-H(\theta_{0})=\Delta Y; \frac{\partial H}{\partial\theta}|_{\theta=\theta_{0}}=J(\theta)$

则可以把关于 $\theta$ 的非线性最小二乘问题化为关于 $\Delta \theta$ 的线性最小二乘问题

$\Delta Y = J(\theta_{0})\Delta \theta+v$

则可以采用线性最小二乘法的方法求得每一次迭代的 $\Delta \theta$

$\Delta \theta=[J^{T}J]^{-1}J^{T}DY$

二、实例

已知数据集数据如下图所示

已知数据点存于data.txt中，其中第一列是x，第二列是y数据。且已知数据点的表达式为下式

$y=\theta_{1} \exp \left(-\frac{1392}{x}\right)+\frac{1}{\theta_{2} \exp \left(-\sqrt{\frac{\theta_{3}}{\theta_{4}+x}}\right)}=H(x,\theta)+0$

其中有未知参数 $\theta_{1},\theta_{2},\theta_{3},\theta_{4}$ .

【问题】通过非线性最小二乘法的高斯牛顿算法对未知参数进行估计，误差保证在 $10^{-3}$ 以内。

【求解】对观测矩阵H求雅可比矩阵

$J(x,\theta))=\left[\begin{array}{llll} \frac{\partial y}{\partial \theta_{1}} & \frac{\partial y}{\partial \theta_{2}} & \frac{\partial y}{\partial \theta_{3}} & \frac{\partial y}{\partial \theta_{4}} \end{array}\right]=\left[\exp \left(-\frac{1392}{x}\right), \frac{-\exp \left(\sqrt{\frac{\theta_{3}}{\theta_{4}+x}}\right)}{\theta_{2}^{2}}, \frac{\exp \left(\sqrt{\frac{\theta_{3}}{\theta_{4}+x}}\right)}{\left.2 \theta_{2} \sqrt{\theta_{3}\left(\theta_{4}+x\right.}\right)}, \frac{-\exp \left(\sqrt{\frac{\theta_{3}}{\theta_{4}+x}}\right)\left(\frac{\theta_{3}}{\theta_{4}+x}\right)^{3 / 2}}{2 \theta_{2} \theta_{3}}\right]$

编写Matlab代码实现高斯牛顿算法对参数迭代求解

%非线性最小二乘法，高斯牛顿法进行迭代求解
clc
clear

theta0 = [1;1;1;1];
data=importdata('data.txt');
length = size(data);

%牛顿迭代法求解
e = 1e-3; thetak=theta0; k=0;
Y = data(:,2);
H = [];       
J = [];

while 1>0
    k = k+1;
    for i=1:length(1)
        H(i) = funy(thetak,data(i,1));
        J(i,:) =  funJ(thetak,data(i,1));
    end
    DY = Y - H';
    Dtheta = inv(J'*J) * J'*DY;
    
    thetak1 = thetak+Dtheta;
    %判断是否跳出迭代
    if Dtheta < e
        break
    end  
    thetak = thetak1;
end

disp('迭代的次数k为:')
disp(k)
disp('近似最优解θ为:')
disp(thetak1)

% tk1=[0.0098;1158.8;11.225;63.892];
yls = [];
for i=1:length
    yls(i) = thetak1(1)*exp(-1392/data(i,1)) +...
             1 / (thetak1(2)*exp(-sqrt(thetak1(3)/(thetak1(4)+data(i,1)))));
end
plot(data(:,1),yls,'*')
hold on 
plot(data(:,1),data(:,2))
hold off
xlabel('x');ylabel('y');


function H=funJ(tk,x)
H = [exp(-1392/x) ...
    -exp(sqrt(tk(3)/(tk(4)+x)))/(tk(2)^2) ...
    exp(sqrt(tk(3)/(tk(4)+x)))  /  ( 2 * tk(2) * sqrt( tk(3) * (tk(4)+x) ) ) ...
    -exp(sqrt(tk(3)/(tk(4)+x))) * (tk(3)/(tk(4)+x))^(3/2) /  (2 * tk(2) * tk(3) ) ];%1*4
end


function y=funy(theta,x)
    y = theta(1)*exp(-1392/x)+1/(theta(2)*exp(-sqrt(theta(3)/(theta(4)+x))));
end

求解结果

命令行窗口：