% MATLAB SCRIPT FILE s322hw2.m  */
%    SOLUTION TO HW 2, ABOUT SIMPLE SIMULATION AND CONFIDENCE INTERVALS
%    DUE SEPTEMBER 5, 1996
note = 'Running MATLAB script file s322hw2.m' 

iquest = 1 ;      %  question number 1,2,3

format compact ;


if iquest == 1 ;
  n = 400 ;
  note = 'Generating Uniforms' 
          %  I LIKE TO LOOK AT THINGS LIKE THIS, TO SEE THAT THE
          %  PROGRAM IS PROGRESSING AS EXPECTED
  u = rand(n) ;
          %  SINCE SQUARE MATRICES ARE SO COMMON, THIS IS 
          %  UNDERSTOOD TO MEAN "n x n" .


  note = 'Start double loop matrix sum' 
  tic ;
  total = 0 ;
          %  INITIALIZE THIS AT 0, THEN ADD AT EACH STEP
  for i = 1:1:n ;
    if rem(i,100) == 0 ;
      progress = 'working on:' 
      i
    end ;

    for j = 1:1:n ;
      total = total + u(i,j) ;
    end ;
  end ;
  outs = toc ;
  outs = ['For double loop, elapsed time (secs) was: ', num2str(outs)] ;
          %  num2str CONVERTS "NUMBERS" TO "STRINGS" OF CHARACTERS

  note = 'output results to file'
  fid = fopen('s322hw2.out','wt') ;
    titlstr1 = 'Output from the MATLAB Script s322hw2.m,   ' ;
    titlstr1 = [titlstr1,date] ;
    titlstr2 = '  Solution to Home Work Assignment 2, due September 5, 1996' ;
  cntbytes = fprintf(fid,'%1s\n',titlstr1) ;
  cntbytes = fprintf(fid,'%1s\n',titlstr2) ;

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


  note = 'Start single loop matrix sum' 
  tic ;
  total = 0 ;
  for i = 1:1:n ;
    if rem(i,100) == 0 ;
      progress = 'working on:' 
      i
    end ;

    total = total + sum(u(:,i)) ;
  end ;
  outs = toc ;
  outs = ['For single loop, elapsed time (secs) was: ', num2str(outs)] ;

  note = 'output results to file'
  cntbytes = fprintf(fid,'%1s\n',outs) ;


  note = 'Start direct matrix sum' 
  tic ;
  total = sum(sum(u)) ;
  outs = toc ;
  outs = [' For direct sum, elapsed time (secs) was: ', num2str(outs)] ;
  
  note = 'output results to file'
  cntbytes = fprintf(fid,'%1s\n',outs) ;

  fclose('all') ;


elseif iquest == 2 ;

  n = 60 ;
  note = 'Generating Uniforms' 
  u = rand(n) ;

  note = 'Start triple loop matrix product' 
  tic ;
  for i = 1:1:n ;
    if rem(i,20) == 0 ;
      progress = 'working on:' 
      i
    end ;

    for j = 1:1:n ;
      product = 0 ;  
          %  INITIALIZE FOR EACH TERM
      for k = 1:1:n ;
        product = product + u(i,k) * u(k,j) ;
          %  BUILD UP SUMMATION FOR EACH TERM
      end ;
    mprod(i,j) = product ;
    end ;
  end ;
  outs = toc ;
  outs = ['For triple loop, elapsed time (secs) was: ', num2str(outs)] ;

  note = 'output results to file'
  fid = fopen('s322hw2.out','at') ;
          %  THIS OPENS FOR "APPENDING", I.E. KEEPING OLD FILE
          %  CONTENTS, AND ADDS NEW STUFF
    headstr1 = 'Problem 2:  ' ;
  cntbytes = fprintf(fid,'\n\n %1s\n',headstr1) ;
  cntbytes = fprintf(fid,'%1s\n',outs) ;


  note = 'Start double loop matrix product' 
  tic ;
  for i = 1:1:n ;
    if rem(i,20) == 0 ;
      progress = 'working on:' 
      i
    end ;

    for j = 1:1:n ;
      mprod(i,j) = u(i,:) * u(:,j) ;  
    end ;
  end ;
  outs = toc ;
  outs = ['For double loop, elapsed time (secs) was: ', num2str(outs)] ;

  note = 'output results to file'
  cntbytes = fprintf(fid,'%1s\n',outs) ;


  note = 'Start direct matrix product' 
  tic ;
  mprod = u * u ;  
  outs = toc ;
  outs = ['For direct prod, elapsed time (secs) was: ', num2str(outs)] ;

  note = 'output results to file'
  cntbytes = fprintf(fid,'%1s\n',outs) ;

  fclose('all') ;


elseif iquest == 3 ;
  nobs = 100 ;
  nds = 10000 ;
          %  NUMBER OF DATA SETS
  ndspb = 100 ;
          %  NUMBER OF DATA SETS PER BLOCK
  nblocks = nds / ndspb ;
          %  NUMBER OF BLOCKS

  note = 'Loop through blocks' 
  tnout = 0 ;
         %  INITIALIZE THIS FOR COUNTING NUMBER OF TIMES OUTSIDE
  for iblock = 1:1:nblocks ;
    if rem(iblock,20) == 0 ;
      progress = 'working on block number:' 
      iblock
    end ;

    data = rand(nobs,ndspb) ;
          %  GENERATE WHOLE BLOCK OF DATA
    vmean = mean(data) ;
          %  VECTOR OF MEANS FOR THIS BLOCK
    vsd = std(data) ;
          %  VECTOR OF S.D.'S FOR THIS BLOCK
    cil = vmean - 1.96 * vsd ./ sqrt(nobs) ;
          %  LEFT ENDPOINTS OF CONFIDENCE INTERVALS
    cir = vmean + 1.96 * vsd ./ sqrt(nobs) ;
          %  RIGHT ENDPOINTS OF CONFIDENCE INTERVALS

    flagl = (.5 < cil) ;
          %  FLAGS (WITH A "1") THOSE ENTRIES IN THE VECTOR CIL,
          %  (AND HENCE THOSE DATA SETS), WHERE THE TRUE MEAN IS
          %  BELOW THE LEFT ENDPOINT
    flagr = (.5 > cir) ;
          %  SIMILAR FOR RIGHT END
    nout = sum(flagl + flagr) ;
          %  SUMMING 0'S AND 1'S IS SAME AS COUNTING NUMBER OF 1'S
          %  I.E. NUMBER OF TIMES MEAN IS OUTSIDE C.I.
    tnout = tnout + nout ;
  end ;


  note = 'output results to file'
  fid = fopen('s322hw2.out','at') ;
    headstr1 = 'Problem 3:  ' ;
  cntbytes = fprintf(fid,'\n\n %1s\n',headstr1) ;

    outs = ['In a simulation run of ', num2str(nds)] ;
    outs = [outs, ' data sets of size ', num2str(nobs)] ;
  cntbytes = fprintf(fid,'%1s\n',outs) ;
    outs = 'The number of times a 95% CI did not contain the ' ;
    outs = [outs, 'mean was: ', num2str(tnout)] ;
  cntbytes = fprintf(fid,'%1s\n',outs) ;
    outs = 'I.e. the empirically observed coverage probability was: ' ;
    outs = [outs, num2str(1 - tnout/nds)] ;
  cntbytes = fprintf(fid,'%1s\n',outs) ;

  fclose('all') ;

end ;
