%MATLAB SCRIPT FILE s322hw4.m  */
%
%    SOLUTION TO HOMEWORK 4 FOR STATISTICS 322
%    POSTED 9/26/96
%    EXPLORES REJECTION METHOD GENERATION OF BETA'S

format compact ;
disp('Running s322hw4.m') ;

iquest = 1 ;    %  1,...,5
                               %  2 DOES OUTPUT FOR 1 & 2
                               %  5 DOES OUTPUT FOR 5 & 6

if iquest == 1 ;        %  THEN ANSWER QUESTION 1: 
                        %        DEVELOP FUNCTION fbeta

  %  NOTE: THIS FILE WAS EDITED SIMULTANEOUSLY WITH fbeta.m
                       
  itest = 1 ;   %  DIFFERENT TESTS TO TRY:  1,...,6

  xtest = linspace(-.1,1.1,1000) ;
          %  VECTOR OF VALUES TO TEST AT

  if itest == 1 ;
    denstr = 'Uniform density, ' ;
    pa = 1 ;
    pb = 1 ;
          %  THE NAME beta IS NOT USED FOR A VARIABLE, AS THAT IS
          %  THE MATLAB COMMAND FOR THE BETA FUNCTION

  elseif itest == 2 ; 
    denstr = 'Triangular density, ' ;
    pa = 2 ;
    pb = 1 ;

  elseif itest == 3 ; 
    denstr = 'Parabolic density, ' ;
    pa = 3 ;
    pb = 1 ;

  elseif itest == 4 ; 
    denstr = 'Pole at zero density, ' ;
    pa = .5 ;
    pb = 1 ;

  elseif itest == 5 ; 
    denstr = 'Asymmetric U-shape density, ' ;
    pa = .8 ;
    pb = .1 ;
          %  WHY IS THE MASS TEST HORRIBLE HERE?
          %  TRY SOME OTHER ALPHAS

  elseif itest == 6 ; 
    denstr = 'Roughly Normal, with mean 1/3, ' ;
    pa = 10 ;
    pb = 20 ;

  end ;
  titstr = [denstr 'alpha = ' num2str(pa) ] ;
  titstr = [titstr ' beta = ' num2str(pb)] ;
  
  density = fbeta(xtest,pa,pb) ;

  %  SEE IF PICTURE LOOKS AS EXPECTED
  disp('Make plot') ;
  plot(xtest,density) ;
    title(titstr) ;
    
  %  DO NUMERICAL INTEGRAL TO SEE IF MASS IS NEARLY ONE
  mass = trapz(xtest,density) ;

    massstr = ['Area under curve ~ ' num2str(mass)] ;
    theight = axis ;
    theight = theight(3) + .9 * (theight(4) - theight(3)) ;
  text(.1,theight,massstr) ;
          %  PUT THE ESTIMATED MASS ON THE PLOT


elseif iquest == 2 ;    %  THEN ANSWER QUESTION 2, 
                        %        WHY WON'T REJECTION WORK FOR ALPHA,
                        %        BETA < 1?

  %  SIMPLY WRITE ANSWER INTO A FILE
  fid = fopen('s322hw4.out','wt') ;
    titlstr1 = 'Output from the MATLAB Script s322hw4.m,   ' ;
    titlstr1 = [titlstr1,date] ;
    titlstr2 = '  Solution to Home Work Assignment 4, ' ;
    titlstr2 = [titlstr2 'due September 26, 1996'] ;
  cntbytes = fprintf(fid,'%1s\n',titlstr1) ;
  cntbytes = fprintf(fid,'%1s\n',titlstr2) ;

    headstr1 = 'Problem 1:  ' ;
    str11 = 'See the MATLAB script fbeta.m, and run s322hw4.m for tests' ;
  cntbytes = fprintf(fid,'\n\n\n    %1s\n',headstr1) ;
  cntbytes = fprintf(fid,'        %1s\n',str11) ;


    headstr2 = 'Problem 2:  ' ;
    str21 = '    This Rectangular box rejection method won''t work' ;
          %  NOTE: USE '' TO GET A ' INTO STRINGS
    str22 = 'for alpha < 1 or beta < 1, because then the density has' ;
    str23 = 'a pole, so the larger region will have infinite area.' ;
  cntbytes = fprintf(fid,'\n\n\n    %1s\n',headstr2) ;
  cntbytes = fprintf(fid,'    %1s\n',str21) ;
  cntbytes = fprintf(fid,'    %1s\n',str22) ;
  cntbytes = fprintf(fid,'    %1s\n',str23) ;

elseif iquest == 3 ;    %  THEN ANSWER QUESTION 3, 
                        %        DEVELOP FUNCTION genbeta

  itest = 1 ;      % 1,2

  if itest == 1 ;    %  VISUAL CHECK BASED ON P-P PLOT
    pa = 1 ;    %  HAVE TRIED 1,2,5,100
    pb = 1 ;    %  HAVE TRIED 1,3,20
    nr = 1000 ;
    nc = 1 ;
    rvout = genbeta(nr,nc,pa,pb) ;

    %  GENERATE COMPONENTS OF P-P PLOT
      pdel = 1 / (nr + 1) ;
    pgrid = linspace(pdel,1-pdel,nr) ;
          %  EQUALLY SPACED GRID OF PROBABILITIES
    ordstat = sort(rvout) ;
          %  ORDER STATISTICS OF DATA

    plot(pgrid,betainc(ordstat,pa,pb)) ;
          %  EMPIRICAL PROBABILITIES, VS. THEORETICAL PROBABILITIES
        titstr = ['P-P plot to test Beta(' num2str(pa)] ;
        titstr = [titstr ',' num2str(pb) ') R. V. s'] ;
      title(titstr) ;
      xlabel('Empirical') ;
      xlabel('Theoretical') ;

  elseif itest == 2 ;   %  CHECK MATRIX VERSION
    pa = 1 ;    %  1,
    pb = 1 ;    %  1,
    nr = 5 ;
    nc = 4 ;
    disp('Here comes a matrix of Beta r.v.s') ;
    testbeta = genbeta(nr,nc,pa,pb) 

  end ;

  fid = fopen('s322hw4.out','at') ;
    headstr3 = 'Problem 3:  ' ;
    str31 = 'See the MATLAB script genbeta.m, and run s322hw4.m for tests' ;
  cntbytes = fprintf(fid,'\n\n\n    %1s\n',headstr3) ;
  cntbytes = fprintf(fid,'        %1s\n',str31) ;


elseif iquest == 4 ;    %  THEN ANSWER QUESTION 4, 
                        %        CONSTRUCT TABLE CHECKING PROBS

  %  WRITE PRELIMINARY STUFF TO OUTPUT FILE
  fid = fopen('s322hw4.out','at') ;
    headstr4 = 'Problem 4:  ' ;
    str41 = '                  x <= .25           x <= .5' ;
    str41 = [str41 '           x <= .75'] ;
    str42 = 'alpha beta    RCI   Th.   LCI    RCI   Th.   LCI' ;
    str42 = [str42 '    RCI   Th.   LCI'] ;
  cntbytes = fprintf(fid,'\n\n\n    %1s\n',headstr4) ;
  cntbytes = fprintf(fid,'%1s\n',str41) ;
  cntbytes = fprintf(fid,'%1s\n',str42) ;

  nr = 1000 ;
  xgrid = [.25, .5, .75] ;
  mxgrid = vec2mat(xgrid,nr) ;
          %  MATRIX VERSION FOR USE BELOW
  mparam = [1,2; 12,4; 20,1; 40,40] ;
          %  MATRIX OF PARAMETER VALUES
  
  %  LOOP THORUGH PAIRS OF PARAMETERS
  for iparam = 1:size(mparam,1) ;
    pa = mparam(iparam,1) ;
    pb = mparam(iparam,2) ;

    data = genbeta(nr,1,pa,pb) ;
    mdata = vec2mat(data,size(xgrid,2)) ;
          %  MATRIX VERSION FOR COMPARISON

    eprobs = mean(mdata < mxgrid) ;
          %  EMPIRICAL PROBABILITIES
    tprobs = betainc(xgrid,pa,pb) ;
          %  THEORETICAL PROBABILITIES, FROM INCOMPLETE BETA FUNCTION

      sds = sqrt(tprobs .* (1 - tprobs)) ;
          %  BINOMIAL STANDARD DEVIATIONS
    cil = eprobs - 1.96 * sds / sqrt(nr) ;
    cir = eprobs + 1.96 * sds / sqrt(nr) ;
          %  LEFT AND RIGHT ENDPOINTS OF CONFIDENCE INTERVALS

    %  WRITE RESULTS TO TABLE
      frmtstr = ' %4.1f %4.1f  %6.3f%6.3f%6.3f' ;
      frmtstr = [frmtstr ' %6.3f%6.3f%6.3f %6.3f%6.3f%6.3f\n'] ;
      vout = [cil; tprobs; cir] ;
      vout = vout(:)' ;
          %  MAKE INTO SINGLE COLUMN VECTOR
    cntbytes = fprintf(fid,frmtstr,[pa pb vout]) ;

  end ;

elseif iquest == 5 ;    %  THEN ANSWER QUESTION 5, 
                        %        CONSTRUCT TABLE OF COMPUTATION TIMES

  nr = 4000 ;
  pa = 1 ;
  pb = 60 ;

  %  WRITE FIRST PART TO OUTPUT FILE
  fid = fopen('s322hw4.out','at') ;
    headstr5 = 'Problem 5:  ' ;
    str51 = 'The 10 computation times were:' ;
  cntbytes = fprintf(fid,'\n\n\n    %1s\n',headstr5) ;
  cntbytes = fprintf(fid,'    %1s\n',str51) ;

  tottime = 0 ;
  for i = 1:10 ;
    tic ;
          %  START STOPWATCH TIMER

    data = genbeta(nr,1,pa,pb) ;

    out = toc ;
          %  STOP STOPWATCH TIMER, AND SAVE IN VECTOR OF TIMES

    tottime = tottime + out ;
          %  KEEP RUNNING TOTAL

    %  OUTPUT TIME FOR THIS RUN
      frmtstr = ['           %12.5f\n'] ;
    cntbytes = fprintf(fid,frmtstr,out) ;

  end ;

  %  OUTPUT TOTAL TIME
    str52 = 'The total computation time was:' ;
  cntbytes = fprintf(fid,'\n    %1s\n',str52) ;
    frmtstr = ['           %12.5f\n'] ;
  cntbytes = fprintf(fid,frmtstr,tottime) ;


  %  OUTPUT ANSWER FOR PROBLEM 6
    headstr6 = 'Problem 6:  ' ;
  cntbytes = fprintf(fid,'\n\n\n    %1s\n',headstr6) ;
    str61 = 'The algorithm is very slow because there is a lot of "waste"' ;
    str62 = 'meaning very many uniforms need to be generated to get a' ;
    str63 = 'sample of beta RVs.' ;
  cntbytes = fprintf(fid,'        %1s\n',str61) ;
  cntbytes = fprintf(fid,'    %1s\n',str62) ;
  cntbytes = fprintf(fid,'    %1s\n',str63) ;


end ;
fclose('all') ;

