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

    function [TS,OK] =pevinit(TS,Xgc,Yrf,W0,univariate,isCoded);
%PEVINIT initializes a two-stage model ready for evaluating pev

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


if nargin<5
	univariate=0;
end
if nargin<6 | ~isCoded
	[Xgc,Yrf,W0,OkGroups,lossdf]= checkdata(TS,Xgc,Yrf,W0);
else
	lossdf=0;
end
if length(TS.covmodel)==0
   TS= covinit(TS,Xgc,Yrf);
end

nf= size(Yrf,2);
nc= length(TS.covmodel);
% degrees of freedom
% nf*ns - p - numDispersionParams
% df= numel(Yrf)-numParams(TS)- nc - length(covmodel(TS.Local)) - lossdf; % local dispersion parameters
df= numel(Yrf)-numParams(TS)-  lossdf; % local dispersion parameters


OK= canMLE(TS) && df>0;

if OK
	
   if univariate
	  % make covariance diagonal
	  sig= unstruct(TS.covmodel);
	  sig= diag(diag(sig));
	  TS.covmodel= covupdate(TS.covmodel,sig);
     	lam= RidgeMatrix(TS);
   else
      lam= [];
   end
   
   % scale data
   [Ts2,Xs,Ys,W0s,S]= mle_scale(TS,Xgc,Yrf,W0);
   
   Wci= choltinv(Ts2.covmodel,W0s);
   % make sure we are doing qr on a full matrix
	

	if nnz(lam)
		Xa= [Wci*Xs;lam];
	else
		Xa= Wci*Xs;
	end
	
	[Xs,S]= xregprecond(Xa);
	R= qr(Xs,0);
    % check rank of R
    rd= abs(diag(R));
    tol= length(Xa)*eps*max(rd);
    if tol==0
        tol= eps;
    end
    OK=  ~(size(Xa,2)>size(R,1) | any(rd<tol));
    if OK
        if issparse(R)
            ri= S/R;
        else
            n= size(R,2);
            ri= S/(triu(R(1:n,:)));
        end

        % store results
        TS= var(TS,ri,[],df);
    end
end