disp(['Running MATLAB script s322eg14.m']) ;
%
%    EXAMPLE 14 FOR STATISTICS 322
%    POSTED 12/5/96
%    EXPLORES FOURIER BASIS NONPARAMETRIC REGRESSION
%    

nobs = 200 ;
sig = .5 ;       %  STANDARD DEVIATION OF ADDITIVE GAUSSIAN NOISE 

left = -1 ;      %  LIMITS FOR BOTH GENERATION AND PLOTTING 
right = 1 ;
top = 1.5 ;
bottom = -1.5 ;
randn('seed',90987234) ;

disp('Generating Data') ;
xgrid = linspace(left,right,nobs)' ;
          %  VECTOR OF x_i POINTS, AS A COLUMN VECTOR
mofx = 1 - abs(2 * xgrid) ;
          %  REGRESSION FUNCTION
y = mofx + sig .* randn(nobs,1) ;
          %  ADD GAUSSIAN NOISE


disp('Save data to a file for later use') ;
diary s322eg14.dat ;
  format long ;
          %  SAVES MORE SIGNIFICANT DIGITS
  [xgrid mofx y]
diary off ;


disp('Fitting Fourier Functions and making plot') ;
vk = [0; 1; 4; 10] ;
          %  VECTOR OF CUTOFF DEGREES k
mhat = [] ;
nk = size(vk,1) ;
for ik = 1:nk ;
  k = vk(ik) ;
  disp(['    Fitting Fourier Series of order ' num2str(k)]) ;

  %  BUILD X MATRIX
  x = sqrt(2) * ones(nobs,1) ;
  for iik = 1:k ;
    x = [x cos(iik * pi * xgrid)] ;
    x = [x sin(iik * pi * xgrid)] ;
  end ;

  %  CALCULATE ESTIMATE
  xpxi = inv(x' * x) ;
  vcoeff = xpxi * x' * y ;
  mhat = [mhat (x * vcoeff)] ;

  vcoeffm = x' * mofx ;
  vcoeffm = xpxi * vcoeffm ;
end ;

subplot(2,1,1) ;
  plot(xgrid,mhat(:,1),'.',...
       xgrid,mhat(:,2),':',... 
       xgrid,mhat(:,3),'--',... 
       xgrid,mhat(:,4),'-.',xgrid,y,'+',xgrid,mofx,'-') ;
%    legend('k = 0','k = 1','k = 4','k = 10','Data','m(x)') ;
          %  DIDN'T LIKE THE LOOKS OF THIS, DECIDED 
          %                      LEGEND BELOW WAS ENOUGH
      title('Data and Fourier Series Est''s') ;
      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),'LineStyle'),'-') ;
          %  THE LINE BELOW DIDN'T WORK HERE, BECAUSE '-' CAME
          %  OUT AS EQUAL TO '--', WHICH GAVE THE WRONG RESULT.
%        if get(vaxchil(ic),'LineStyle') == '-' ;
          %  CHANGE THE LINE THICKNESS FOR THE CURVE WITH
          %  THE SOLID LINE STYLE
          %  NOTE THIS STRUCTURE AUTOMATICALLY CHOOSES THE CORRECT
          %  CURVE, EVEN IF THE PROGRAM IS CHANGED LATER
          set(vaxchil(ic),'LineWidth',2.0) ;
          disp('Changed Line Width') ;
          %  WRITE THIS TO SCREEN, TO MAKE SURE WENT OK
        end ;
        if strcmp(get(vaxchil(ic),'LineStyle'),'.') ;
          set(vaxchil(ic),'MarkerSize',3) ;
          disp('Changed Marker Size') ;
        end ;
      end ;


disp('Do ctrl-brk to stop and create diary file') ;
pause
%  DID CTRL-BRK HERE TO STOP AND RUN COMMANDS THAT WERE SAVED IN
%  DIARY FILE s322eg14.dia.  THOSE WERE USED TO LEARN HOW TO 
%  CHANGE LINE THICKNESSES, ETC. ABOVE.

subplot(2,1,2) ;
  plot(xgrid,mhat(:,1),'.',...
       xgrid,mhat(:,2),':',... 
       xgrid,mhat(:,3),'--',... 
       xgrid,mhat(:,4),'-.',xgrid,mofx,'-') ;
    legend('k = 0','k = 1','k = 4','k = 10','m(x)',1) ;
    title('Fourier Series Est''s') ;
    axis([-1 1 -1.1 1.1]) ;
    vaxchil = get(gca,'Children') ;
    for ic = 1:length(vaxchil) ;     %  LOOP THROUGH CURVES
      if strcmp(get(vaxchil(ic),'LineStyle'),'-') ;
        set(vaxchil(ic),'LineWidth',2.0) ;
        disp('Changed Line Width') ;
      end ;
      if strcmp(get(vaxchil(ic),'LineStyle'),'.') ;
        set(vaxchil(ic),'MarkerSize',3) ;
        disp('Changed Marker Size') ;
      end ;
    end ;

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




disp('Any key for bar graph') ;
pause ;

subplot(2,3,1) ;
  bar(vcoeff.^2) ;
  title('Individ. Powers, Data') ;
  axis([0 22 0 .8]) ;

subplot(2,3,2) ;
  bar(vcoeffm.^2) ;
  title('Individ. Powers, Signal') ;
  axis([0 22 0 .8]) ;

subplot(2,3,3) ;
  bar((vcoeff - vcoeffm).^2) ;
  title('Individ. Powers, Noise') ;
  axis([0 22 0 .01]) ;

subplot(2,3,4) ;
  bar(6 + log10(vcoeff.^2)) ;
  title('6 + log10 Powers, Data') ;
  axis([0 22 0 6]) ;

subplot(2,3,5) ;
  bar(6 + log10(vcoeffm.^2)) ;
  title('6 + log 10 Powers, Signal') ;
  axis([0 22 0 6]) ;

subplot(2,3,6) ;
  bar(6 + log10((vcoeff - vcoeffm).^2)) ;
  title('6 + log 10 Powers, Noise') ;
  axis([0 22 0 6]) ;

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

