disp(['Running MATLAB script s322hw8.m']) ;
%
%    SOLUTION TO HOMEWORK 8 FOR STATISTICS 322
%    POSTED 10/31/96
%    EXPLORES LEGENDRE ORTHONORMAL BASIS
format compact ;

ipart = 1 ;    %  1 - TEST CREATION OF MATRIX OF BASIS FUNCTIONS
               %  2 - MAKE PLOT OF BASIS FUNCTIONS
               %  3 - TEST COEFFICIENTS
               %  4 - MAKE PLOT OF PROJECTIONS

ia = 1 ;       %  1,2,3  (WHEN ipart = 3,4) 


va = [1; 4; 10] ;
left = - 1 ;
right = 1 ;

if ipart == 1 | ipart == 3 ;    %  THEN DO TESTING 
  nxgrid = 10000 ;   
          %  INTERESTING TO TRY BOTH 100 AND 1000 HERE
  nmax = 4 ;
elseif ipart == 2 | ipart == 4 ;    %  THEN MAKE PLOTS
  nxgrid = 50 ;   
  nmax = 4 ;
end ; 

xgrid = linspace(left,right,nxgrid)' ;
xinc = xgrid(2) - xgrid(1) ;

disp('Create Matrix of Polynomials') ;
if nmax == 0 ;    %  THEN ONLY NEED CONSTANT VECTOR
  mp = ones(nxgrid,1) ;
else ;            %  THEN BUILD UP BIGGER MATRIX
  mp = [ones(nxgrid,1) xgrid] ;
  n = 2 ;
  while n <= nmax ;   %  CONTINUE TO ADD COLUMNS
    p = (2 * n - 1) .* xgrid .* mp(:,n) ;
    p = p - (n - 1) * mp(:,n-1) ;
    p = p / n ;
    mp = [mp p] ;
    n = n + 1 ;
  end ;
end ;

mpsi = vec2mat(sqrt(.5:1:(nmax+.5)),nxgrid) .* mp ;

if ipart == 1 ;
  disp('Doing Tests of Basis Functions') ;
  format long ;

    psi0 = sqrt(0+.5) * ones(nxgrid,1) ;
    e0 = max(abs(mpsi(:,1) - psi0)) ;
  disp(['If 0th basis function is right, this is 0: ' num2str(e0)]) ;

    psi1 = sqrt(1+.5) * xgrid ;
    e1 = max(abs(mpsi(:,2) - psi1)) ;
  disp(['If 1st basis function is right, this is 0: ' num2str(e1)]) ;

    psi2 = sqrt(2+.5) * .5 * (3 * xgrid.^2 - 1) ;
    e2 = max(abs(mpsi(:,3) - psi2)) ;
  disp(['If 2nd basis function is right, this is 0: ' num2str(e2)]) ;

    psi3 = sqrt(3+.5) * .5 * (5 * xgrid.^3 - 3 * xgrid) ;
    e3 = max(abs(mpsi(:,4) - psi3)) ;
  disp(['If 3rd basis function is right, this is 0: ' num2str(e3) 10 10]) ;

  format short ;
  disp(['Matrix of inner products (should look like identity '...
                     'matrix) is:']) ;
  matinp = (xinc * (mpsi' * mpsi))


elseif ipart == 2 ;
  disp('Making the Plot for Part 2') ;

  titstr = 'Statistics 322 HW #8, Part 5' ;

  plot(xgrid,mpsi(:,1),'o',xgrid,mpsi(:,2),'+',...
       xgrid,mpsi(:,3),':',xgrid,mpsi(:,4),'--',...
       xgrid,mpsi(:,5),'-') ;
    title(titstr) ;
    legend('n=0','n=1','n=2','n=3','n=4',0) ;

  print -dps hw8a.ps ;

elseif ipart >= 3 ;
  disp('Get Coefficients') ;
  a = va(ia) ;
  f = exp(a * xgrid) ;

  %  TO CREATE APPROX'S, FIRST GET COEFFICIENTS
  vij = 2 * sinh(a) / a ;
          %  FIRST Ij,  I0
  vcoeff = sqrt(1/2) * vij ;
          %  FIRST (I.E. 0TH) COEFFICIENT
  mapprox = mpsi(:,1) * vcoeff ;
          %  FIRST APPROXIMATION
  for n = 1:nmax ;
    np1 = n + 1 ;
          %  BECAUSE VECTORS ARE INDEXED STARTING WITH 1, NOT 0
    ij = 2 * sinh(a) / a ;
    ij = ij - 4 * n * cosh(a) / a.^2 ;
    ij = ij + 2 * n * (2 * n - 1) * vij(n) / a.^2 ;
    vij = [vij; ij] ;

    vj = (0:1:n)' ;
          %  VECTOR OF J'S IN SUM FOR COEFFICIENT
          %  NEXT START BUILDUP OF THE COEFFICIENT
    coeff = gamma(n + 1) ./ (gamma(vj + 1) .* gamma(n - vj + 1)) ;
          %  RECALL gamma(n + 1) = n!
    coeff = coeff .* (-1).^(n - vj) ;
    coeff = coeff .* vij ;
    coeff = sum(coeff) ;
    coeff = sqrt(n + .5) * (-a).^n * coeff / (2^n * gamma(n + 1)) ;
    vcoeff = [vcoeff; coeff] ;

    mapprox = [mapprox (mpsi(:,1:np1) * vcoeff)] ;
          %  MULTIPLY EACH COLUMN BY APPROPRIATE COEFFICIENT    

  end ;

  if ipart == 3 ;
    disp('Doing Tests of Coefficients') ;
    format long ;

    disp([10 10 'Numerically Integrated Ij']) ;
    niij = (xinc * (vec2mat(xgrid',nmax+1).^(2*vec2mat(vj,nxgrid)) * f)) ;
    disp('         Exact          Approximated       |Difference|') ;
    [vij niij abs(vij - niij)]

    disp([10 10 'Numerically Integrated Coefficients']) ;
    nicoeff = (xinc * (mpsi' * f)) ;
    disp('         Exact          Approximated       |Difference|') ;
    [vcoeff nicoeff abs(vcoeff - nicoeff)]

  elseif ipart == 4 ;
    disp('Making the Plot for Part 4') ;

    titstr = 'Statistics 322, HW #8, Part 8' ;
    titstr = [titstr ', a = ' num2str(a)] ;

    normf = sinh(2 * a) / a ;
    fpower = cumsum(vcoeff.^2) / normf ;

    plot(xgrid,f,'-',xgrid,mapprox(:,1),'.',...
         xgrid,mapprox(:,2),':',xgrid,mapprox(:,3),'--',...
         xgrid,mapprox(:,4),'o',xgrid,mapprox(:,5),'*');
    legend('f',['n=0, %power=' num2str(fpower(1))],...
               ['n=1, %power=' num2str(fpower(2))],...
               ['n=2, %power=' num2str(fpower(3))],...
               ['n=3, %power=' num2str(fpower(4))],...
               ['n=4, %power=' num2str(fpower(5))],0) ;
    title(titstr) ;

    eval(['print -dps hw8b' num2str(a) '.ps']) ;

  end ;

end ;

