disp('MATLAB script file s322hw6.m') ;
%
%    SOLUTION TO HOMEWORK 6 FOR STATISTICS 322
%    POSTED 10/10/96
%    EXPLORES BOOTSTRAP CONFIDENCE INTERVALS

iquest = 2 ;         %  QUESTION NUMBER  2,3

if iquest == 2 ;     %  THEN ANSWER QUESTION 2

  vn = [10, 25, 100] ;
  nn = size(vn,2) ;
            %  VECTOR OF SAMPLE SIZES
  nsim = 1000 ;
            %  NUMBER OF SIMULATIONS
  nboot = 499 ;
            %  NUMBER OF BOOTSTRAP SAMPLES
  
  rand('seed',92369742) ;
  randn('seed',30485703) ;
            %  USE SEEDED GENERATION, FOR EXACTLY REPRODUCIBLE RESULTS


  mcp1 = zeros(nn,3)  ; 
  mcp2 = zeros(nn,3)  ; 
  mcp3 = zeros(nn,3)  ; 
            %  MATRICES FOR STORING FINAL COVERAGE PROBABILITIES

  for in = 1:nn ;
            %  LOOP THROUGH SAMPLE SIZES
    n = vn(in) ;

    for idist = 1:3 ;

      %  FIRST GENERATE MATRIX OF DATA, WITH EACH COLUMN THE DATA FOR
      %  ONE SIMULATION.
      if idist == 1 ;
        datstr = 'N(0,1) data' ;
        mdata = randn(n,nsim) ;
        mu = 0 ;
      elseif idist == 2 ;
        datstr = 'Expon(1) data' ;
        mdata = -log(rand(n,nsim)) ;
        mu = 1 ;
      elseif idist == 3 ;
        datstr = 'Chi2(1) data' ;
        mdata = randn(n,nsim).^2 ;
        mu = 1 ;
      end ;

      %  GET VECTORS OF SUMMARY STATISTICS 
      vxbar = mean(mdata)' ;
            %  COLUMN VECTOR IS NICEST FOR WORK BELOW 
      vs = std(mdata)' ;
            %  1/(n-1) VERSION
      vsdhat = sqrt((n - 1) / n) * vs ;
            %  1/n VERSION

      %  DO ORDINARY CENTRAL LIMIT THEOREM CI'S AS MATRIX OPERATION
      ztable = 1.282 ;
      vci1 = vxbar - ztable * vsdhat / sqrt(n) ;

      %  DO STUDENT'S t TABLE CI'S AS MATRIX OPERATION
      if n == 10 ;
        ttable = 1.383 ;
      elseif n == 25 ;
        ttable = 1.318 ;
      elseif n == 100 ;
        ttable = 1.294 ;
      end ;
      vci2 = vxbar - ttable * vs / sqrt(n) ;

      %  NOW DO LOOP FOR BOOTSTRAP PART
      vqstar1 = [] ;
      vqstar2 = [] ;
      for isim = 1:nsim ;
        disp([datstr ', n = ' num2str(n) ', sim # ' num2str(isim)]) ;

        vdata = mdata(:,isim) ;
        randind = ceil(n * rand(n,nboot)) ;
            %  RANDOM INDICES FROM 1,...,N
        mbdata = vdata(randind) ;
            %  COLUMN VECTOR OF BOOTSTRAP DATA
        mbdata = reshape(mbdata,n,nboot) ;
            %  MATRIX OF BOOTSTRAP DATA
      
        xbboot1 = sqrt(n) * (mean(mbdata)' - vxbar(isim)) ;
            %  NORMALIZED VERSION OF BOOTSTRAP XBAR
          sdboot = sqrt((n - 1) / n) * std(mbdata)' ;
        xbboot2 = xbboot1 ./ sdboot ;
            %  STUDENTIZED VERSION OF BOOTSTRAP XBAR

        xbboot1 = sort(xbboot1) ;
        xbboot2 = sort(xbboot2) ;
            %  SORT, SO INDICES EASILY GIVE QUANTILES

        vqstar1 = [vqstar1; xbboot1(round(.90 * (nboot + 1)))] ;
        vqstar2 = [vqstar2; xbboot2(round(.90 * (nboot + 1)))] ;
            %  GET .90 QUANTILES, AND FORM INTO COLUMN VECTORS

      end ;

      vci3 = vxbar - vqstar1 / sqrt(n) ;
      vci4 = vxbar - vqstar2 .* vsdhat / sqrt(n) ;
            %  LOWER 90% BOOTSTRAP CONFIDENCE BOUNDS

      %  CALCULATE COVERAGE PROBABILITIES
      mcp1(in,idist) = mean(vci1 < mu) ;
      mcp2(in,idist) = mean(vci2 < mu) ;
      mcp3(in,idist) = mean(vci3 < mu) ;
      mcp4(in,idist) = mean(vci4 < mu) ;
            %  MEAN OF 0'S AND 1'S GIVES EMPIRICAL FREQUENCY OF 1'S

    end ;

  end ;


  %  GET MATRICES OF STANDARD ERRORS
  mse1 = sqrt(mcp1 .* (1 - mcp1) / nsim) ;
  mse2 = sqrt(mcp2 .* (1 - mcp2) / nsim) ;
  mse3 = sqrt(mcp3 .* (1 - mcp3) / nsim) ;
  mse4 = sqrt(mcp4 .* (1 - mcp4) / nsim) ;


  %  WRITE ANSWERS INTO A FILE
  fid = fopen('s322hw6.out','wt') ;
    titlstr1 = 'Output from the MATLAB Script s322hw6.m,   ' ;
    titlstr1 = [titlstr1,date] ;
    titlstr2 = '  Solution to Home Work Assignment 6, ' ;
    titlstr2 = [titlstr2 'due October 10, 1996'] ;
  cntbytes = fprintf(fid,'%1s\n',titlstr1) ;
  cntbytes = fprintf(fid,'%1s\n',titlstr2) ;

    headstr2 = '  Problem 2:  ' ;
  cntbytes = fprintf(fid,'\n\n\n%1s\n',headstr2) ;

    str21 = '    For the N(0,1) distribution:' ;
    str22 = '            Stand. C.L.T.     Student''s t     Naive Boot. ' ;
    str22 = [str22, '     Stud. Boot.'] ;
    str23 = '             C.P.  (s.e.)     C.P.  (s.e.)    C.P.  (s.e.)' ;
    str23 = [str23, '     C.P.  (s.e.)'] ;
  cntbytes = fprintf(fid,'\n%1s\n',str21) ;
  cntbytes = fprintf(fid,'%1s\n',str22) ;
  cntbytes = fprintf(fid,'%1s\n',str23) ;
    fstr = ' n = %3.0f    %5.3f (%5.3f)    %5.3f (%5.3f)    ' ;
    fstr = [fstr '%5.3f (%5.3f)    %5.3f (%5.3f)\n'] ;
    outm =[vn; mcp1(:,1)'; mse1(:,1)'; mcp2(:,1)'; mse2(:,1)'; ...
                    mcp3(:,1)'; mse3(:,1)'; mcp4(:,1)'; mse4(:,1)'] ;
  cntbytes = fprintf(fid,fstr,outm) ;

    str21 = '    For the Exp(1) distribution:' ;
  cntbytes = fprintf(fid,'\n%1s\n',str21) ;
  cntbytes = fprintf(fid,'%1s\n',str22) ;
  cntbytes = fprintf(fid,'%1s\n',str23) ;
    outm =[vn; mcp1(:,2)'; mse1(:,2)'; mcp2(:,2)'; mse2(:,2)'; ...
                    mcp3(:,2)'; mse3(:,2)'; mcp4(:,2)'; mse4(:,2)'] ;
  cntbytes = fprintf(fid,fstr,outm) ;

    str21 = '    For the Chi2(1) distribution:' ;
  cntbytes = fprintf(fid,'\n%1s\n',str21) ;
  cntbytes = fprintf(fid,'%1s\n',str22) ;
  cntbytes = fprintf(fid,'%1s\n',str23) ;
    outm =[vn; mcp1(:,3)'; mse1(:,3)'; mcp2(:,3)'; mse2(:,3)'; ...
                    mcp3(:,3)'; mse3(:,3)'; mcp4(:,3)'; mse4(:,3)'] ;
  cntbytes = fprintf(fid,fstr,outm) ;

  fclose('all') ;


elseif iquest == 3 ;     %  THEN ANSWER QUESTION 3
  disp('Working on Question 3') ;

  n = 10 ;
  nboot = 499 ;
  rand('seed',6836424) ;
  randn('seed',1909374) ;

  for idist = 1:3 ;

    %  FIRST GENERATE VECTOR OF DATA
    if idist == 1 ;
      datstr = 'N(0,1)' ;
      vdata = randn(n,1) ;
      mu = 0 ;
      sd = 1 ;
    elseif idist == 2 ;
      datstr = 'Expon(1)' ;
      vdata = -log(rand(n,1)) ;
      mu = 1 ;
      sd = 1 ;
    elseif idist == 3 ;
      datstr = 'Chi2(1)' ;
      vdata = randn(n,1).^2 ;
      mu = 1 ;
      sd = sqrt(2) ;
    end ;
    xbar = mean(vdata)' ;

    %  NOW DO BOOTSTRAP PART (COPIED FROM ABOVE)
    randind = ceil(n * rand(n,nboot)) ;
          %  RANDOM INDICES FROM 1,...,N
    mbdata = vdata(randind) ;
          %  COLUMN VECTOR OF BOOTSTRAP DATA
    mbdata = reshape(mbdata,n,nboot) ;
          %  MATRIX OF BOOTSTRAP DATA
    xbboot1 = sqrt(n) * (mean(mbdata)' - xbar) ;
          %  NORMALIZED VERSION OF BOOTSTRAP XBAR
      sdboot = sqrt((n - 1) / n) * std(mbdata)' ;
    xbboot2 = xbboot1 ./ sdboot ;
          %  STUDENTIZED VERSION OF BOOTSTRAP XBAR
    xbboot1 = xbboot1 ./ sd ;
          %  PUT ON SCALE OF n(0,1) 

      delb = .4 ;
    vbinc = -4:delb:4 ;
          %  VECTOR OF HISTOGRAM BIN CENTERS
    gpdf = (1 / sqrt(2 * pi)) * exp(-vbinc.^2) ;
          %  GAUSSIAN PROBABILITY DENSITY FUNCTION
    gpdf = nboot * delb * gpdf ;
          %  NORMALIZED TO SCALE OF COUNTS

    subplot(3,3,3*(idist - 1) + 1) ;
      hist(xbboot1,vbinc) ;
        hold on ;
          plot(vbinc,gpdf,'b') ;
        hold off ;
        axis([-4.1,4.1,0,100]) ;
        title([datstr ' Direct']) ;

    subplot(3,3,3*(idist - 1) + 2) ;
      hist(xbboot2,vbinc) ;
        hold on ;
          plot(vbinc,gpdf,'b') ;
        hold off ;
        axis([-4.1,4.1,0,100]) ;
        title([datstr ' Student']) ;

    subplot(3,3,3*(idist - 1) + 3) ;
      hist(sdboot,0:delb/4:2) ;
        axis([0,2,0,200]) ;
        title([datstr ' Sig hat']) ;

  end ;

  print -dps s322hw6.ps

end ; 
