www.gusucode.com > mbcmodels 工具箱 matlab 源码程序 > mbcmodels/@xreglinear/qrdecomp.m

    function [Q,R,OK,df,varargout]= qrdecomp(m, x, Wc, J)
%QRDECOMP QR decomposition for least squares
%
% [Q,R,OK,df, Ri]= qrdecomp(m,x,Wc,J);
% Inputs
%   m  model
%   x observation matrix
%   Wc optional weights
%   J  jacobian
% Outputs
%   Q,R economised qr decomposition of QR=J , R square upper triangular
%       note this is adjusted for lambda and type of least squares (e.g. ridge, rols)
%   OK rank check on R (J)
%   df degrees of freedom
%   Ri is defined as factorisation of  Ri*Ri'= (R'*R)\J'*J/(R'*R)

%  Copyright 2000-2007 The MathWorks, Inc. and Ford Global Technologies, Inc.

if nargin < 4
    % compute the Jacobian if we are not passing it in
    J = CalcJacob(m,x);
    if nargin < 3
        Wc = [];
    end
end

if isempty(J)
    df= size(J,1);
    Q=zeros(df,0);
    R= [];
    OK=df>0;
    varargout= cell(1,nargout-4);
else
    switch m.qr
        case 'ols'
            [Q,R,OK,df,varargout{1:nargout-4}]= qr_ols(m,x,Wc,J);
        case 'ridge'
            [Q,R,OK,df,varargout{1:nargout-4}]= qr_ridge(m,x,Wc,J);
        case 'rols'
            [Q,R,OK,df,varargout{1:nargout-4}]= qr_rols(m,x,Wc,J);
        otherwise
            [Q,R,OK,df,varargout{1:nargout-4}]= feval(['qr_',m.qr],m,x,Wc,J);
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% QR decomposition for Ordinary Least Squares
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Q,R,OK,df, Ri]= qr_ols(m,x,Wc,J)

if ~isempty(Wc)
    % weight the Jacobian
    J = Wc*J;
end

[Q,R,OK]= xregqr(J);

df = size(Q,1)-size(R,1);

if nargout > 4
    if OK
        Ri = inv(R);
    else
        Ri = [];
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% QR decomposition for ROLS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Q,R,OK,df,Ri]= qr_rols(m,x,Wc,J)

if ~isempty(Wc)
    % weight the Jacobian
    J = Wc*J;
end

if numel(m.lambda)==1
    lam= m.lambda;
else% local rols
    termsin = Terms(m);
    if any(termsin)
        lam= m.lambda(termsin);
        lam= lam(:);
    else
        lam = 0;
    end
end


n = size(J,2);

% use normalised qr decomposition
[Q,R,OK]= xregqr(J);
if OK
    dB= diag(R);
    Db = sparse(1:n,1:n,dB,n,n);
    W= Q*Db;
    B= Db\R;

    % W is not orthonormal
    dw= (sum(W.*W,1))';
    dwlam= sqrt(lam+dw);

    % normalise to get Q,R equivalents so b= R\(Q'*y) and H= sum(Q.*Q,2)
    % note Q'*Q ~= I
    Ddwlam = sparse(1:n,1:n,1./dwlam,n,n);
    Q= W*Ddwlam;
    R= Ddwlam\B;
else
    df=0;
    Ri=[];
    return
end


df = size(J,1) - sum(dw./(dw+lam));% N - gamma

if nargout > 4
    if OK
        % inverse calculation for PEV and std(b) = Ri*Ri'
        D = sparse(1:n,1:n,sqrt(dw./(lam+dw)),n,n);
        Ri= R\D;
    else
        Ri = [];
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% QR decomposition for RIDGE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Q,R,OK,df,Ri]= qr_ridge(m,x,Wc,J)

if ~isempty(Wc)
    % weight the Jacobian
    J = Wc*J;
end

if numel(m.lambda)==1
    sqrtlam= sqrt(m.lambda)*eye(size(J,2));
else % local ridge
    termsin = Terms(m);
    if any(termsin)
        sqrtlam= diag(sqrt(m.lambda(termsin)));
    else
        sqrtlam = eye(size(J,2));
    end
end

OK = size(J,1)>=size(J,2);
if OK

    % augmented J
    Ja = [J; sqrtlam];

    % apply Q, R to the augmented matrix
    [Q,R,OK]= xregqr(Ja);
    Q= Q(1:size(J,1),:);

    dq = sum(Q.*Q,2);
    df = size(J,1) - sum(dq);

    if nargout > 4
        if OK
            Ri = R\Q';
        else
            Ri = [];
        end
    end
else
    Q=[];
    R=[];
    df=0;
    Ri=[];
end