基于最小二乘算法的参数估计（matlab程序和测试）

最小二乘法是参数估计中最常见的一种算法，在时间序列ARMA、曲线拟合、神经网络、最优化等领域有着非常广泛的应用。

这里我分享一下我写的一个最小二乘估计一个MATLAB函数，里面实现的算法有三种：递推最小二乘算法、带遗忘因子的最小二乘算法、偏差补偿的最小二乘算法。下面先给出一些测试结果，同时也是对用法的一个说明。

测试使用的时间序列函数为： z(k+2) = 1.3*z(k+1) - 0.5*z(k) + 0.8*u(k+1) + 0.5*u(k) +v(k)

首先生成测试序列z(k)和u(k)：

%Z(k+2)=1.3*Z(k+1)-0.5*Z(k)+0.8*u(k+1)+0.5*u(k)+v(k)clearclcx = [0 1 0 1 1 0 1 1 1];n = 400;                u = [];                  for i=1:n    temp = xor(x(4),x(9));    u(i) = x(9);    for j = 9:-1:2        x(j) = x(j-1);    end    x(1) = temp;endv = randn(1,398);z = zeros(400,1);z(1) = -1;z(2) = 0;for i=3:400z(i) = 1.5*z(i-1)-0.7*z(i-2)+u(i-1)+0.5*u(i-2)+v(i-2);end

生成的z(k)序列如下图所示

图1 生成的测试序列z(k)

测试一：递推最小二乘参数估计

调用函数来进行参数估计

>> a = RLS(2,2,z,u,'RLS',[],[]);

得到的结果是

图2 使用递推最小二乘参数估计得到的收敛过程

得到参数估计值

待估计的参数值为：-1.54330.7231161.047250.370914

结果正确。

测试二： 带遗忘因子的最小二乘参数估计

调用函数

>> a = RLS(2,2,z,u,'fRLS',0.99,2);

结果如下

图3 带遗忘因子的参数估计

图3可以明显看到带遗忘因子的参数估计的收敛曲线波动较大。但是这种算法对参数是时变的时间序列参数估计很有效。

图4 带遗忘因子的参数估计的估计方差收敛曲线

得到的参数

待估计的参数值为：-1.546030.7250451.15610.369839

测试三：偏差补偿最小二乘递推算法

a = RLS(2,2,z,u,'oRLS',[],1)

结果：

误差图不画了，占空间。

得到的参数估计值

待估计的参数值为：-1.561510.7464631.234270.360575

以上就是一个简单的测试过程，也说明了函数的调用方法。

最后给出函数

function theta = RLS(M,N,z,u,method,beta,display)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Parameter estimation using recursive least squre method% 使用递推最小二乘法估计ARMA模型参数,返回估计参数值% ==========================参数说明===================================% 输出参数(return parameters)：%      theta [vec] -  返回ARMA参数估计% 输入参数(input parameters):%         M  [num] -  AR阶数%         N  [num] -  MA阶数%         z  [vec] -  输入数据，用于估计参数%         u  [vec] -  MA的噪声序列，实际上是很难得到准确的序列的%    method  [cha] -  使用的最小二乘法类型%            - 'RLS'：递推最小二乘法，The recursive least square method%            - 'fRLS'：带遗忘因子的最小二乘递推算法%            - 'oRLS'：偏差补偿最小二乘递推算法%      beta  [num] -  遗忘因子%   display  [num] -  画出估计参数参数变化曲线和误差变化曲线%            - '0' ：不画图%            - '1' ：只画出参数估计变化曲线%            - '2' : 同时画出参数变化曲线和误差变化曲线% =============================syntax==================================% 在使用最小二乘估计ARMA参数之前，应该准备好观测序列z,MA噪声序列u% 如果没有u,将自动视为AR参数估计。% example1 ： a = RLS(2,2,z,u,'RLS',[],[])%             使用递推最小二乘参数估计，AR和MA阶数分别为2,2,后两项参数为空 % example2  : a = RLS(2,2,z,u,'fRLS',0.99,2)%             使用带遗忘因子的最小二乘递推算法，遗忘因子设为0.99%             同时画出估计参数的收敛过程和估计误差的收敛过程% example3 ： a = RLS(2,2,z,u,'oRLS',[],1)%             使用补偿最小二乘估计算法% ======================================================================% Written by zhangwenyu,2013,12,12% If you have any suggestions, mail me at [e-mail]%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%n = length(z);len = M+N;% check argumentsif isempty(method)    method = 'RLS';endif isempty(beta)    beta = 0.98;endif isempty(display)    display = 1;endif M <= 1    error('输入参数中ar的阶数M必须大于1');endif M < N    error('设定M的值大于N，请修改！');endif isempty(u)    Mk = 0.5*randn(1,n+1)./std(randn(1,n+1));     %白噪声    warning('MA噪声序列自动设定为白噪声');else    Mk = u;end% A Mk (colored noise) generate process % x = [0 1 0 1 1 0 1 1 1];   % for i=1:n+len-1%     temp=xor(x(4),x(9));%     Mk(i)=x(9);%     for j=9:-1:2%         x(j)=x(j-1);%     end%     x(1)=temp;% end% Initialization settingsP = 2*eye(len);Pstore = zeros(len,n-1);for k=1:len    Pstore(k,1) = P(k,k);endtheta = zeros(len,n-1);theta(:,1) = 2;h = zeros(len,1);% choose which method to useswitch method % RLS method        case 'RLS'    for i = M+1:n        for p=1:M            h(p) = -z(i-p);        end        for q=1:N            h(M+q) = Mk(i-q);        end        K = P*h*inv(h'*P*h+1);        theta(:,i-1) = theta(:,i-2)+K*(z(i)-h'*theta(:,i-2));        P = (eye(len)-K*h')*P;        for j=1:len             Pstore(j,i-1) = P(j,j);        end    end% RLS with forgetting factor% useful in time-varying parameter estimation       case 'fRLS'    for i = M+1:n       for p=1:M           h(p) = -z(i-p);       end       for q=1:N           h(M+q) = Mk(i-q);       end       K = P*h*inv(h'*P*h+1);       theta(:,i-1) = theta(:,i-2)+K*(z(i)-h'*theta(:,i-2));       P = (eye(len)-K*h')*P/beta;       for j=1:len            Pstore(j,i-1) = P(j,j);       end      end% RLS with bias compensation method    case 'oRLS'        thetaO = zeros(len,n-1);        thetaO(:,1) = 2;        D = [eye(M) zeros(M,N);zeros(N,len)];        J = 0;    for i = M+1:n       for p=1:M           h(p) = -z(i-p);       end       for q=1:N           h(M+q) = Mk(i-q);       end       J=J+(z(i-1)-h'*theta(:,i-1))^2/(1+h'*P*h);       K = P*h*inv(h'*P*h+1);       theta(:,i-1) = theta(:,i-2)+K*(z(i)-h'*theta(:,i-2));       P = (eye(len)-K*h')*P/beta;       for j=1:len            Pstore(j,i-1) = P(j,j);       end       es=J/((i-1)*(1+(thetaO(:,i-2))'*D*theta(:,i-1)));       thetaO(:,i-1)=theta(:,i-1)+(i-1)*es*P*D*thetaO(:,i-2);    endend% set display options    i=1:n-1;      switch display        case 0            fprintf('不画图\n');        case 1            figure;            for k=1:len                plot(i,theta(k,:));                hold on;            end            title('待估参数的收敛过程');        case 2            figure;            for k=1:len                plot(i,theta(k,:));                hold on;            end            title('待估参数的收敛过程');            figure;            for k=1:len                plot(i,Pstore(k,:));                hold on;            end                       title('估计方差的收敛过程');    end% print estimated parameters in command windowfprintf('待估计的参数值为：\n');for k=1:len    fprintf('%g\n',theta(k,end));endend