www.gusucode.com > mbcmodels 工具箱 matlab 源码程序 > mbcmodels/@xregrbf/wigglecenters.m
function [om,OK]=wigglecenters(m) %WIGGLECENTERS Place centers where model is wiggly % % WIGGLECENTERS is a novel algorithm to move the centers to where y is % wiggly. % Copyright 2000-2014 The MathWorks, Inc. and Ford Global Technologies, Inc. om= contextimplementation(xregoptmgr,m,@i_wigglecenters,[],'WiggleCenters',@wigglecenters); om = setAltMgrs(om,{@wigglecenters,@rols,@rederr}); % give a choice of center selection algorithms % fit parameters om= AddOption(om,'MaxNCenters','min(nObs/4,25)','evalstr', 'Number of centers');% percentage of data taken as centers om= AddOption(om,'PercentCandidates','(min(nObs,200)/nObs)*100','evalstr', 'Percentage of data to be candidate centers');% number of candidate centers,take min(nObs,200) om= AddOption(om,'Tolerance',0.0001,{'numeric',[eps 1]}, [],false);% stopping criterion om= AddOption(om,'cost',0,'numeric',[],false); OK = 1; function [m,cost,OK]=i_wigglecenters(m,om,x0,x,y,varargin) set(m,'qr','rols'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % set maxncenters nObs = size(x,1); maxCStr= get(om,'MaxNCenters'); maxncenters = calcMaxNCenters(m,maxCStr,nObs); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% nf = get(m,'nfactors'); %%%%%%%%%%%adjust the number of candidate centers in pool PcCandStr = get(om,'PercentCandidates'); PercentCandidates = calcPercentCand(m,PcCandStr,nObs); approxncand = max(1,round((PercentCandidates/100)*nObs)); %approximate number of candidate points to try if nObs > approxncand candcenters = xregrbfcentersel(approxncand,x); else % take all data points as candidate centers candcenters = x; end m.centers = candcenters; ncand = size(candcenters,1); %%%%%%%%%%% % Get the parameters stored in m width = m.width;%width of the radial basis function if any(width < eps) || size(width,1)> 1 %covers when on entry we have a width per center [m,OK] = defaultwidth(m,x);%set the default width end lambda=get(m,'lambda'); if length(lambda)>1 % only use one lambda lambda = lambda(end); set(m,'lambda',lambda); end % Set up the 'full' regression matrix, Phik with all the candidate centers % At the kth step Phik will store columns corresponding to the centers that % have yet to be selected. % The columns will be orthogonal to the (k-1) columns corresponding centers % already chosen Phik = x2fx(m,x); k=1; radius = sqrt(nf)/maxncenters; distances = xregdistancepoints(candcenters,x); neighbors = (distances<radius); numneigh = sum((neighbors),1);% number of neighbors within radius centorder = 1:ncand; d = y; % will store the residual y2 = y'*y; B = zeros(maxncenters,ncand); %initialise the upper triangular matrix B g= zeros(1,min(ncand,maxncenters)); while k<= min(ncand,maxncenters), %loop on k B(k,:)= 0; B(k,k)= 1; drep = repmat(d,1,ncand-k+1);% yneigh = zeros(nObs,ncand-k+1); yvarneigh = zeros(nObs,ncand-k+1); %step 1 %compute the regularised weight (gk) tmp = sum(Phik.*Phik,1) + lambda; gk = (d'*Phik)./tmp; % yvalues at the neighbors of the candidate centers yneigh(neighbors) = drep(neighbors); ave = sum(yneigh,1)./numneigh;% average of the y values of neighbours ave = repmat(ave,nObs,1); yvar = yneigh-ave; yvarneigh(neighbors) = yvar(neighbors); wiggle = sum(yvarneigh.^2,1)'; [mxwiggle,kindex] = max(wiggle); % choose the candidate with the largest wiggle index = kindex + k-1; % interchanging the selected column of Phik with the 1st column of Phik Phik(:,[1 kindex]) = Phik(:,[kindex 1]); neighbors(:,[1 kindex]) = neighbors(:,[kindex 1]); numneigh([1 kindex]) = numneigh([kindex 1]); % interchange (up to the (k-1)th row) the selected column of B with the kth column of B B(1:k-1,[k index]) = B(1:k-1,[index k]); % interchange the order of the centers centorder([k index])= centorder([index k]); % get the kth regularised weight g(k) = gk(kindex); % name the first column of Phik (corresponds to the kth center) wk = Phik(:,1); %kill the first column to remove this center from consideration Phik(:,1) = []; neighbors(:,1) = []; numneigh(1) = []; %Step 3 %Gram-Schmidt orthogonalisatiom to make the %remaining columns of Phi orthogonal to wk ww = wk'*wk; B(k,k+1:ncand) = wk'*Phik/ww; Phik = Phik - (wk/ww)*(wk'*Phik); % (this formula appears in regularisation in the selection...) % Step 4 % compute the remaining residual d = d - g(k)*wk; k = k+1; end nchosen = k-1;%number of centers selected %save the centers chosen. %note centers are in the order that they were chosen m.centers = candcenters((centorder(1:nchosen)),:); B = B(1:nchosen,1:nchosen);% the final upper triangular matrix weights = B\g(1:nchosen)';%get the unregularised weights m = update(m,weights);%load the weights into the xreglinear parent cost = log10GCV(m,x,y);% setFitOpt(m,'cost',cost); OK=1;