disp(['Running MATLAB script s322eg15.m']) ;
%
%    EXAMPLE 15 FOR STATISTICS 322
%    POSTED 12/7/96
%    EXPLORES REGRESSION SPLINE FITS
%    

vdeg = [0; 1; 2; 3; 4; 8] ;
          % COLUMN VECTOR OF DEGREES TO SHOW

nx = 50 ;
left = 0 ;
right = 1 ;
bottom = -1.5 ;
top = 1.5 ;
xgrid = linspace(left,right,nx)' ;
y = sin(2 * pi * xgrid) ;

vknot = (0:(1/3):1)' ;
          %  COLUMN VECTOR OF KNOTS

for ivdeg = 1:6 ;
  deg = vdeg(ivdeg) ;
  disp(['   Working on degree ' num2str(deg)]) ;

  %  CONSTRUCT THE DESIGN MATRIX
  %  POLYNOMIAL PART FIRST
  x = ones(nx,1) ;
  if deg > 0 ;
    for ideg = 1:deg ;
      x = [x (xgrid.^ideg)] ;
          %  FILL OUT MATRIX WITH COLUMNS CORRESPONDING TO POWERS
    end ;
  end ;

  %  KNOT PART SECOND
  mramp = [] ;
  for iknot = 1:(length(vknot)-2) ;
    if deg == 0 ;
      ramp = xgrid - vknot((iknot+1),1) ;
      ramp = (ramp > 0) ;
    else ;
      ramp = xgrid - vknot((iknot+1),1) ;
      flag = (ramp > 0) ;
      ramp = ramp .^ deg ;
      ramp = flag .* ramp ;
    end ;
    x = [x ramp] ;
    mramp = [mramp ramp] ;
  end ;

  %  CALCULATE LEAST SQUARES
  coeff = inv(x' * x) * x' * y ;
  mfit = x * coeff ;

  %  BUILD MATRIX OF POLY PIECES TO PLOT
  poly = x(:,1:(deg+1)) * coeff(1:(deg+1)) ;
  mpoly = poly ;
  for ipiece = 2:(length(vknot) - 1) ;
    poly = poly + coeff(deg+ipiece) * (xgrid - vknot(ipiece)).^deg ;
    mpoly = [mpoly poly] ;
  end ;


  subplot(3,2,ivdeg) ;
    plot(xgrid,y,'+w',xgrid,mfit,'-y',xgrid,mpoly,'.r') ;
      axis([left right bottom top]) ;
      title('Regression Spline Fits') ;
      text(.75,1.2,['Deg = ' num2str(deg)]) ;
      hold on ;
        plot(vec2mat(vknot(2:3)',2),vec2mat([top; bottom],2),':b') ;
          %  ADD ON VERTICAL LINES SHOWING KNOTS
      hold off ;
      vaxchil = get(gca,'Children') ;
          %  VECTOR OF OBJECT HANDLES FOR CHILDREN OF AXES, WHICH
          %  ARE THE CURVES.
      for ic = 1:length(vaxchil) ;     %  LOOP THROUGH CURVES
        if strcmp(get(vaxchil(ic),'Type'),'line') ;
          if strcmp(get(vaxchil(ic),'LineStyle'),'-') ;
          %  CHANGE THE LINE THICKNESS FOR THE CURVE WITH
          %  THE SOLID LINE STYLE
            set(vaxchil(ic),'LineWidth',1.0) ;
            disp('          Changed Line Width') ;
          end ;
          if strcmp(get(vaxchil(ic),'LineStyle'),'.') ;
            set(vaxchil(ic),'MarkerSize',3) ;
            disp('          Changed Marker Size') ;
          end ;
        end ;
      end ;

end ;

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