www.gusucode.com > mbcdesign 工具箱 matlab 源码程序 > mbcdesign/@des_multimod/vdelete.m

    function [des,psi]=vdelete(des,initpsi,p,DO_DESIGNTYPE)
%DES_LINEARMOD/VDELETE   V-optimal deletion
%   [D,PSI]=VDELETE(D,INITPSI,P) deletes P lines from the design D
%   using v-optimality.  A new design object and the new
%   v-opimality criteria, PSI, are returned.

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



% Created 12/6/2000

if nargin<4
   DO_DESIGNTYPE=1;
end

[des,K] = sumxtx(des);
fp = freepoints(des);
fs = factorsettings(des);

mmdl=model(des);
% set up matrices for each model
mdls = get(mmdl,'models');
wts = get(mmdl,'weights');
nmdls = get(mmdl,'nmodels');
X=cell(1,nmdls);
Ai=cell(1,nmdls);
for n=1:nmdls
   % Initial Ai must be for all points.  It is then iteratively
   % updated and hence is always correctly derived from all points.
   X{n}=CalcJacob(mdls{n},fs);
   
   [~,R]= qr(X{n},0);
   ri= R\eye(size(R));
   
   Ai{n}= ri*ri';
   
   % X used for deletion point search is taken from the non-fixed
   % points in the design
   X{n}=X{n}(fp,:);
end

if DO_DESIGNTYPE
   % need to remember start point
   [TP,INFO]= DesignType(des);
end

psinew= initpsi;

for m=1:p
   psiupd = zeros(length(fp)-m+1,1);
   for n=1:nmdls
      % build up div as weighted sum
      % quick way
      B= X{n}*Ai{n};
      div=(1-sum(B.*X{n},2));
      
      Bi=wts(n).*sum( B.*(B*K{n}),2 );
      % check div isn't 0!!
      div=max(div,Bi*eps*2);
      % form output vector
      psiupd = psiupd + Bi./div;
   end
   
   [delpsi,i] = min(psiupd);   
   psinew = psinew+delpsi;
   
   if p>1
      for n=1:nmdls
         % Update Ai
         % Need to recalc the ith parts for each model
         Bi = X{n}(i,:)*Ai{n};
         divi = (1-sum(Bi.*X{n}(i,:),2));
         
         Ai{n}= mx_r1update(Ai{n},Bi/sqrt(divi),0);
         
         X{n}(i,:)=[];
      end
   end
   des=delete(des,'indexed',i,'changeable');
end
if DO_DESIGNTYPE
   % update design type
   des=DesignType(des,TP,INFO);         % reset object to initial setting
   des=UpdateDesignType(des,'v');       % update type setting
end

psi=psinew;
return