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