ridgecg.m

偏最小二乘算法在MATLAB中的实现

function [b,thetamin,cumpress] = ridgecg(xblock,yblock,tmaxm,divs,split,gui,action);
%RIDGECG Ridge regression by cross validation for use with MODLGUI

%  Copyright
%  Eigenvector Technologies
%  Graphics modified November 1993
%  Modified May 1995
global press cumpress theta b tmaxm divs split thetamin incr
handl = get(gui,'UserData');
[mx,nx] = size(xblock);
[my,ny] = size(yblock);
if strcmp(action,'crval')
  press = zeros(split,divs+1);
  if mx ~= my
    error('Number of samples must be the same in both blocks')
  end
  incr  = tmaxm/divs;
  for i = 1:split
    ind = zeros(mx,1);
    count = 0;
    for j = 1:mx
      test = round((j+i-1)/split) - ((j+i-1)/split);
      if test == 0
        ind(j,1) = 1;
        count = count + 1;
      end
    end
    [a,b] = sort(ind);
    newx  = xblock(b,:);
    newy  = yblock(b,:);
    calx  = newx(1:mx-count,:);
    testx = newx(mx-count+1:mx,:);
    caly  = newy(1:mx-count,:);
    testy = newy(mx-count+1:mx,:);
    [m,n] = size(xblock);
    dfs   = m - n - 1;
    b     = zeros(n,divs+1);
    b(:,1) = calx\caly;
    ridi  = diag(diag(calx'*calx));
    for j = 1:divs
      b(:,j+1) = inv(calx'*calx + ridi*j*incr)*calx'*caly;
    end
    for j = 1:divs+1
      ypred      = testx*b(:,j);
      press(i,j) = sum(sum((ypred-testy).^2));
    end
    figure(handl(1,1))
    plot(0:incr:tmaxm,press(i,:))
    s = sprintf('PRESS for Test Set Number %g out of %g',i,split);
    title(s);
    xlabel('Value of theta');
    ylabel('PRESS');
    drawnow
  end
  pause(2)
  cumpress = sum(press);
  figure(handl(1,1))
  plot(0:incr:tmaxm,cumpress)
  title('Cumulative PRESS as a Function of Theta')
  xlabel('Value of Theta')
  ylabel('PRESS')
  drawnow
  [a,minlv] = min(cumpress);
  thetamin  = (minlv-1)*incr;
  s = sprintf('Minimum Cumulative PRESS is at theta = %g',thetamin);
  figure(gui)
  set(handl(2,4),'String',s)
  set(handl(4,12),'Visible','On')
  set(handl(3,3),'UserData',1);
elseif strcmp(action,'coefs')
  ridi = diag(diag(xblock'*xblock));
  b    = inv(xblock'*xblock + ridi*thetamin)*xblock'*yblock;
  figure(handl(1,1))
  plot(b), hold on, plot(b,'o'), plot(0:nx,zeros(1,nx+1),'-g')
  s = sprintf('Final Regression Vector for Theta = %g',thetamin);
  title(s)
  xlabel('Variable Number')
  ylabel('Regression Coefficient')
  hold off
  figure(gui)
  set(handl(3,3),'UserData',4);
  set(handl(4,12),'Visible','On');
  set(handl(2,4),'String',['Click RESUME to Continue']);
end