disp(['Running MATLAB script s322hw10.m']) ;
%
%    SOLUTION TO HOMEWORK 10 FOR STATISTICS 322
%    POSTED 12/3/96
%    EXPLORES B - SPLINES IN NONPARAMETRIC REGRESSION
%    

format compact ;

load s322eg14.in ;

disp('Check data read in:') ;
data = s322eg14 

xgrid = data(:,1) ;
nx = length(xgrid) ;
m = data(:,2) ;
y = data(:,3) ;
nobs = length(y) ;

left = -1 ;
right = 1 ;
top = 2 ;
bottom = -2 ;

deg = 3 ;

vnknot = [2; 3; 5; 9; 17] ;

disp('  Constructing Design Matrix') ;
sx = ones(nx,1) ;
ideg = 1 ;
while ideg <= deg ;    %  CONSTRUCT THE POLYNOMIAL PART
  sx = [sx (xgrid .^ ideg)] ;
  ideg = ideg + 1 ;
end ;

mbspl = [] ;    %  will store Bsplines
vknot2 = [] ;   %  will store 2nd knots
for inknot = 1:length(vnknot) ;     %  loop through knot numbers
  nknot = vnknot(inknot) ;
  disp(['    Working on ' num2str(nknot) ' knots']) ;

  vknot = linspace(left, right, nknot) ;
  vknot2 = [vknot2 vknot(2)] ;

  x = sx ;    %  START IWTH POLY PARTS ALREADY CONSTRUCTED ABOVE
  for iknot = 1:(length(vknot) - 2) ;
    if deg < 1 ;
      ramp = xgrid - vknot(iknot+1) ;
      ramp = ramp > 0 ;
    else ;
      ramp = xgrid - vknot(iknot+1) ;
      flag = ramp > 0 ;
      ramp = (flag .* ramp).^deg ;
    end ;
    x = [x ramp] ;
  end ;

  disp('        Calculating Least Squares') ;
  coeff = inv(x' * x) * x' * y ;
  mbspl = [mbspl (x*coeff)] ;

end ;


disp('Start Graphics') ;
plot(xgrid,m,'-',xgrid,y,'+',xgrid,mbspl(:,1),'-.',...
        xgrid,mbspl(:,2),'-',xgrid,mbspl(:,3),'--',...
        xgrid,mbspl(:,4),':',xgrid,mbspl(:,5),'.') ;
    titstr = ['Solution to Homework # 10, (printed ' date ')'] ;
  title(titstr) ;
  axis([left right bottom top]) ;
  legend('m(x)','data','2 knots','3 knots','5 knots','9 knots','17 knots') ;

  hold on ;     %  to add vertical lines to show knot locations
   
    for inknot = 2:length(vnknot) ;     %  loop through knot numbers
      nknot = vnknot(inknot) ;
      vknot = linspace(left, right, nknot) ;
      vknot = vknot(2:(nknot-1)) ;      %  interior knots
      if inknot == 2 ;
        plot([vknot; vknot],vec2mat([top; bottom],length(vknot)),'-') ;
      elseif inknot == 3 ;
        plot([vknot; vknot],vec2mat([top; bottom],length(vknot)),'--') ;
      elseif inknot == 4 ;
        plot([vknot; vknot],vec2mat([top; bottom],length(vknot)),':') ;
      elseif inknot == 5 ;
          ndot = 100 ;     %  number of dots, each line
        plot(vec2mat(vknot,ndot),...
               vec2mat(linspace(bottom,top,ndot)',length(vknot)),'.') ;
      end ;
    end ;
  hold off ;

  %  Now text for vertical lines
  text(vknot2(2)+.01,bottom + .95 * (top - bottom),'3 knots') ;
  text(vknot2(3)+.01,bottom + .90 * (top - bottom),'5 knots') ;
  text(vknot2(4)+.01,bottom + .85 * (top - bottom),'9 knots') ;
  text(vknot2(5)+.01,bottom + .80 * (top - bottom),'17 knots') ;

  %  now change line widths
  vaxchil = get(gca,'Children') ; 
          %  vector of handles to children of axes
  naxchil = length(vaxchil) ;
  for iac = 1:naxchil ;
          %  loop through objects on figure
    objh = vaxchil(iac) ;
    if iac == naxchil ;   %  then have target curve
      set(objh,'Linewidth',3) ;
    elseif iac < naxchil - 1 ;  % then look at all others, besides raw data
      if strcmp(get(objh,'Type'),'line') ;
                                  %  then are looking a line, so apply below
        if length(get(objh,'Xdata')) == nobs ; 
                                    %  then have a spline piece
          set(objh,'Linewidth',1) ;
          set(objh,'MarkerSize',5) ;
        elseif length(get(objh,'Xdata')) == ndot ;
                                    %  then have a dotted vertical line
          set(objh,'MarkerSize',2) ;
        end ;
      end ;
    end ;
  end ;

  %  Now add comment text
  text(.3,1.8,'2 & 3 knots not flexible enough') ;
  text(.3,1.6,'17 knots is way too variable') ;
  text(.3,1.4,'5 & 9 knots are rather similar,') ;
  text(.3,1.2,'  personally prefer 5 (more') ;
  text(.3,1.0,'               parsimonious)') ;

  orient landscape ;
eval(['print -dps \matlab\steve\ps\s322hw10.ps']) ;
