%MATLAB SCRIPT FILE s322hw5.m  */
%
%    SOLUTION TO HOMEWORK 5 FOR STATISTICS 322
%    POSTED 10/3/96
%    EXPLORES RANDOM WALKS AND GRAPHICS

disp('Running s322hw5.m') ;

iquest = 1 ;    %  1,2

autop = 0 ;     %  1 FOR PRINTING AS POSTSCRIPT FILES

if iquest == 1 ;        %  THEN ANSWER QUESTION 1: 
                        %        MAKE PLOTS FOR SIMPLE RANDOM WALKS

  rand('seed',92369742) ;
          %  SEEDED GENERATION IS USEFUL, SO RERUNS AREN'T
          %  DISTRACTINGLY DIFFERENT.

  npath = 4 ;
  nmax = 100 ;
  mbern = rand(nmax,npath) ;
          %  LET PATHS RUN DOWN COLUMNS, BECAUSE THIS IS HOW cumsum WORKS.
  mbern = 1 - 2 * (mbern < .5) ;
          %  RANDOM +1'S AND -1'S 
  mpath = cumsum(mbern) ;
          %  GET THE "PARTIAL SUMS" FOR EACH COLUMN

  subplot(2,1,1) ;
  plot(1:nmax,mpath(:,1),'-') ;
  hold on ;
    plot(1:nmax,mpath(:,2),'--') ;
    plot(1:nmax,mpath(:,3),':') ;
    plot(1:nmax,mpath(:,4),'-.') ;
      legend('Path 1','Path 2','Path 3','Path 4') ;
          %  NOTE HOW MATLAB DOES A GOOD JOB OF AUTOMATICALLY
          %  POSITIONING THIS
  hold off ;
    title('Stat 322 HW5, Part 1a') ;
 

  subplot(2,1,2) ;
  stairs(mpath(:,1)) ;
          %  PLOT FIRST PATH
    hold on ;
    for ipath = 2:npath ;
      stairs(mpath(:,ipath)) ;
          %  PLOT NEXT PATH
    end ;
    hold off ;
      title('Stat 322 HW5, Part 1b') ;
 

  if autop == 1 ;
    print s322hw51.ps -dps ;
  end ;     

elseif iquest == 2 ;    %  THEN ANSWER QUESTION 2: 
                        %        MAKE PLOTS FOR 2-d RANDOM WALKS

  rand('seed',69876954) ;
  randn('seed',23245983) ;
  nmax = 10000 ;


  subplot(3,2,1) ;   %  2-D BERNOULLI
    mbern = (rand(nmax,2) < .5) ;
          %  RANDOM 0'S AND 1'S
    mbern = 1 - 2 * (mbern < .5) ;
          %  RANDOM +1'S AND -1'S 
    mpath = cumsum(mbern) ;
          %  GET THE "PARTIAL SUMS"
  plot(mpath(:,1),mpath(:,2)) ;
    title('Stat 322 HW5, Part 2a') ;
    axis([-100 25 -20 160]) ;
          %  DIDN'T LIKE WHAT THE DEFAULT DID, SO CHANGED TO THIS


  subplot(3,2,2) ;   %  3-D BERNOULLI
    vtime = 1:nmax ;
  plot3(mpath(:,1),mpath(:,2),vtime) ;
          %  THIS SHOWS "TIME" ON THE THIRD AXIS, BUT IS A LITTLE
          %  HARD TO SEE, AS IT IS NOT ALWAYS SO CLEAR WHERE THE
          %  PATH IS LYING.  HENCE, ADD VERTICAL LINES, EVERY 1000-TH
          %  TIME POINT.
    hold on ;
      vind = find(rem(vtime,500) == 0) ;
          %  VECTOR OF INDICES FOR EVERY 1000-TH POINT
      xl = mpath(vind,1)' ;
      yl = mpath(vind,2)' ;
      plot3([xl; xl],[yl; yl],[vtime(vind); ones(size(vind))]) ;
          %  THIS ADDS VERTICAL LINES FROM THE "BOTTOM" UP TO THE
          %  PATH.  LENGTH OF THE LINE SHOWS HOW LONG THE TIME HAS
          %  BEEN.
      plot3([0;0],[0;0],[1;nmax],':') ;
          %  PUT A VERTICAL AXIS AT (0,0)
      text(0,0,nmax-1,'time 0') ;
          %  LABEL AXIS TO DISTINGUISH IT FROM OTHERS 
    hold off ;
      view([15000,-2500,10000]) ;
          %  SINCE THE DEFAULT VIEW LETS SOME OF THE VERTICAL LINES
          %  OVERLAP.  THIS VIEW IS CHOSEN BY SOME TRIAL AND ERROR.
      title('Stat 322 HW5, Part 2b') ;
      xlabel('x') ;
      ylabel('y') ;
      axis([-100 25 -20 160 0 nmax]) ;


  subplot(3,2,3) ;   %  2-D GAUSSIAN
    mpath = cumsum(randn(nmax,2)) ;
  plot(mpath(:,1),mpath(:,2)) ;
    title('Stat 322 HW5, Part 2c') ;
      axis([-120 30 -300 30]) ;


  subplot(3,2,4) ;   %  2-D CAUCHY
    mcauch = randn(nmax,2) ./ randn(nmax,2) ;
          %  RATIOS OF NORMALS IS CAUCHY
    mpath = cumsum(mcauch) ; 
  plot(mpath(:,1),mpath(:,2)) ;
    title('Stat 322 HW5, Part 2d') ;
      xmin = 8700 ;
      xmax = 11600 ;
      ymin = -2500 ;
      ymax = 1100 ;
    hold on ;
      plot([xmin xmax xmax xmin xmin], [ymin ymin ymax ymax ymin],'c') ;
            %  THIS DRAWS A "BOX", CORNER BY CORNER, STARTING AT THE
            %  LOWER RIGHT, WHICH WILL SHOW WHERE THE ZOOMED PIECE 
            %  COMES FROM
    hold off ;


  subplot(3,2,5) ;   %  ZOOM 1
  plot(mpath(:,1),mpath(:,2)) ;
    title('Stat 322 HW5, Part 2e') ;
    axis([xmin xmax ymin ymax])
      xmin = 9860 ;
      xmax = 10040 ;
      ymin = -130 ;
      ymax = 150 ;
    hold on ;
      plot([xmin xmax xmax xmin xmin], [ymin ymin ymax ymax ymin],'c') ;
            %  THIS DRAWS A "BOX", CORNER BY CORNER, STARTING AT THE
            %  LOWER RIGHT, WHICH WILL SHOW WHERE THE ZOOMED PIECE 
            %  COMES FROM
    hold off ;


  subplot(3,2,6) ;   %  ZOOM 2
  plot(mpath(:,1),mpath(:,2)) ;
    title('Stat 322 HW5, Part 2f') ;
    axis([xmin xmax ymin ymax])

  if autop == 1 ;
    print s322hw52.ps -dps ;
  end ;     

end ;

