disp(['Running MATLAB script s322eg9.m']) ;
%
%    EXAMPLE 9 FOR STATISTICS 322
%    POSTED 10/8/96
%    EXPLORES MONTE CARLO MARKOV CHAIN

format compact ;

ipart = 2 ;      %  1,2

if ipart == 1 ;    %  THEN DO BASIC CALCULATIONS

  %  WRITE ANSWERS INTO A FILE
  fid = fopen('s322eg9.out','wt') ;
    titlstr1 = 'Output from the MATLAB Script s322eg9.m,   ' ;
    titlstr1 = [titlstr1,date] ;
    titlstr2 = '    Investigations of Monte Carlo Markov Chains, ' ;
  cntbytes = fprintf(fid,'%1s\n',titlstr1) ;
  cntbytes = fprintf(fid,'%1s\n\n\n',titlstr2) ;


  for ijoint = 1:4 ;

    if ijoint == 1 ;
      f = [.10, .14, .10; .10, .10, .13; .13, .10, .10] ;
           %  JOINT DISTRIBUTION MATRIX WITH i,j ENTRY = P[X1=i AND X2=j]
           %  CHOSEN SO THAT ALL ENTRIES ARE RELATIVELY LARGE, 
           %  I.E. GOOD CHANCE OF GOING FROM ANY STATE TO ANY OTHER.
      fstr = 'Good Circulation'
    elseif ijoint == 2 ;
      f = [.30, .02, .02; .02, .28, .02; .02, .02, .30] ;
           %  JOINT DISTRIBUTION MATRIX WITH i,j ENTRY = P[X1=i AND X2=j]
           %  CHOSEN SO THAT EACH STATE IS "STICKY", I.E. SMALL PROB.
           %  OF GOING FROM ANY STATE TO ANY OTHER
      fstr = 'Sticky States' ;
    elseif ijoint == 3 ;
      f = [.01, .16, .16; .16, .02, .16; .16, .16, .01] ;
           %  JOINT DISTRIBUTION MATRIX WITH i,j ENTRY = P[X1=i AND X2=j]
           %  CHOSEN FOR "FAST CIRCULATION", I.E. VERY LARGE PROB.
           %  OF LEAVING CURRENT STATE, BUT GOOD "MIXING" IN VERY RANDOM
           %  CHOICE OF WHERE TO GO NEXT
      fstr = 'Jump Away Fast' ;
    elseif ijoint == 4 ;
      f = [.01, .32, .01; .01, .01, .31; .31, .01, .01] ;
           %  JOINT DISTRIBUTION MATRIX WITH i,j ENTRY = P[X1=i AND X2=j]
           %  CHOSEN FOR "NEARLY CYCLIC", I.E. VERY LARGE PROB.
           %  OF LEAVING CURRENT STATE, TO NEXT ONE
      fstr = 'Nearly Cyclic' ;
    end ;


%    sum(sum(f)) 
          %  THIS WAS A QUICK CHECK THAT THE PROBABILITIES SUM TO 1,
          %  AND THUS WERE E.G. TYPED IN CORRECTLY

    f1 = sum(f')' ;   %  COL VECTOR OF MARGINAL PROBABILITIES FOR X1
    f2 = sum(f) ;     %  ROW VECTOR OF MARGINAL PROBABILITIES FOR X2

    f12 = f ./ vec2mat(f2,3) ;
          %  TRANSITION MATRIX FOR DISTRIBUTION OF X1 GIVEN X2,
          %  WITH i,j ENTRY = P[X1=i | X2=j]
          %      RECALL THE NEW FUNCTION vec2mat 


    f21 = (f ./ vec2mat(f1,3))' ;
          %  TRANSITION MATRIX FOR DISTRIBUTION OF X2 GIVEN X1,
          %  WITH i,j ENTRY = P[X2=i | X1=j]
          %  FINAL TRANSPOSE CHANGES THE ORDER OF i,j

%    disp('Joint probability matrix is:') ;
%    f
%    disp('Marginal Probs for X1 are:') ;
%    f1
%    disp('Marginal Probs for X2 are:') ;
%    f2
%    disp('Transition matrix for X1 given X2 is:') ;
%    f12
%    disp('Transition matrix for X2 given X1 is:') ;
%    f21
%    disp('Any key to continue') ;
%    pause ;
          %  THESE LINES WERE ALSO USED FOR CHECKING, NOW COMMENTED OUT
  
    %  WRITE TO FILE
      outstr = ['    For the ' fstr ' Joint Distribution:'] ;
    cntbytes = fprintf(fid,'\n%1s\n',outstr) ;
      outstr = ['        Joint Probability Matrix is:'] ;
    cntbytes = fprintf(fid,'\n%1s\n',outstr) ;
      fstr = '%5.4f    %5.4f    %5.4f\n' ;
    cntbytes = fprintf(fid,fstr,f) ;

      outstr = ['        Marginal Probs for X1 are:'] ;
    cntbytes = fprintf(fid,'\n%1s\n',outstr) ;
      fstr = '%5.4f    %5.4f    %5.4f\n' ;
    cntbytes = fprintf(fid,fstr,f1) ;


    a = f12 * f21 ;
          %  TRANSITION MATRIX FOR X1 MARKOV CHAIN

%    sum(a) 
          %  THIS WAS USED TO CHECK THAT a REALLY WAS A TRANSITION MATRIX

%    a * f1 
%    f1
          %  USED THIS TO CHECK THAT f1 IS INDEED PRESERVED BY a


    %  NOW ACTUALLY WATCH CONVERGENCE OF PROB'S
      outstr = ['        Sequence of MCMC Probs for X1 are:'] ;
    cntbytes = fprintf(fid,'\n%1s\n',outstr) ;
    format long ;
    vf = [.000001; 0; .999999] ;
    currentvf = vf'  
          %  THIS PART WRITES TO SCREEN
      fstr = 'step = %2.0f   %15.12f   %15.12f   %15.12f\n' ;
    cntbytes = fprintf(fid,fstr,[0 currentvf]) ;
          %  THIS PART WRITES TO FILE
    if ijoint == 1 ;
      mvf1 = currentvf ;
    elseif ijoint == 2 ;
      mvf2 = currentvf ;
    elseif ijoint == 3 ;
      mvf3 = currentvf ;
    elseif ijoint == 4 ;
      mvf4 = currentvf ;
    end ;
          %  BUILD UP THIS WAY, SINCE NO 3D ARRAYS ALLOWED

    for istep = 1:15 ;
      vf = a * vf ;
      currentvf = vf'  
          %  THIS PART WRITES TO SCREEN
        fstr = 'step = %2.0f   %15.12f   %15.12f   %15.12f\n' ;
      cntbytes = fprintf(fid,fstr,[istep currentvf]) ;
          %  THIS PART WRITES TO FILE
      if ijoint == 1 ;
        mvf1 = [mvf1; currentvf] ;
      elseif ijoint == 2 ;
        mvf2 = [mvf2; currentvf] ;
      elseif ijoint == 3 ;
        mvf3 = [mvf3; currentvf] ;
      elseif ijoint == 4 ;
        mvf4 = [mvf4; currentvf] ;
      end ;
          %  BUILD UP THIS WAY, SINCE NO 3D ARRAYS ALLOWED
    end ;

    disp('Any key to continue') ;
    pause ;



    [meigvec,meigval] = eig(a) ;
          %  MEIGVEC IS A MATRIX WITH THE EIGENVECTORS AS COLUMNS
          %  MEIGVAL IS A DIAGONAL MATRIX, WITH THE CORRESPONDING
          %       EIGENVALUES AS ENTRIES

%    meigvec
%    meigval
          %  USED THESE TO SEE RESULT, WAS SURPRISED TO SEE THAT
          %  THE COLUMN OF MEIGVEC CORRESPONDING TO THE 1 EIGENVALUE
          %  WAS not THE SAME AS f1

%    col1 = meigvec(:,1) ;
%    sum(col1.^2) 
          %  THIS SHOWS WHY, EIGENVECTORS ARE NORMALIZED (RECALL ONLY
          %  "DIRECTION", NOT "LENGTH", IS IMPORTANT) TO HAVE EUCLIDEAN
          %  LENGTH = 1.

%    col1 / sum(col1) 
%    f1
          %  USED THESE LINES TO CHECK THAT GET BACK TO f1, BY
          %  NORMALIZING SO THAT SUM OF ENTRIES IS ONE.

%    meigvec' * meigvec 
          %  THIS WAS A SURPRISE FOR ME.  I WAS EXPECTING THIS TO BE
          %  BE THE IDENTITY MATRIX, BUT THAT ONLY HOLDS FOR a 
          %  SYMMETRIC (NOTE NOT FAR FOR SECOND EXAMPLE)

%    max(max(abs(a * meigvec - meigvec * meigval))) 
          %  THESE SHOULD BE 0, BY "EIGENVALUE REPRESENTATION"

    veigval = diag(meigval)' ;
          %  USEFUL COMMAND FOR WORKING WITH DIAGONAL MATRICES
          %  NOTE THIS WORKS "BOTH WAYS", TURNING VECTORS INTO
          %  DIAGONAL MATRICES AND VICE VERSA (DEPENDING ON THE INPUT).

    veigval  
          %  THIS PART WROTE TO SCREEN
      outstr = ['        Eigenvalues for X1 --> X1 transition are:'] ;
    cntbytes = fprintf(fid,'\n%1s\n',outstr) ;
      fstr = '    %8.4f   %8.4f   %8.4f\n' ;
    cntbytes = fprintf(fid,fstr,veigval) ;
          %  THIS PART WROTE TO FILE


  %  NOW LOOK AT EIGENVALUES OF X1 --> X2 TRANSITION MATRIX
    [mevec12,meval12] = eig(f12) ;
    veval12 = diag(meval12)' 
      outstr = ['        Eigenvalues for X1 --> X2 transition are:'] ;
    cntbytes = fprintf(fid,'\n%1s\n',outstr) ;
      fstr = '    %8.4f   %8.4f   %8.4f\n' ;
    cntbytes = fprintf(fid,fstr,veval12) ;


  %  NOW LOOK AT EIGENVALUES OF X2 --> X1 TRANSITION MATRIX
    [mevec21,meval21] = eig(f21) ;
    veval21 = diag(meval21)' 
      outstr = ['        Eigenvalues for X1 --> X2 transition are:'] ;
    cntbytes = fprintf(fid,'\n%1s\n',outstr) ;
      fstr = '    %8.4f   %8.4f   %8.4f\n\f' ;
            %  NOTE /f GIVES A PAGE FEED
    cntbytes = fprintf(fid,fstr,veval21) ;


    disp('Any key to continue') ;
    pause ;
  end ;

  fclose('all') ;

elseif ipart == 2 ;    %  THEN DO GRAPHICS

  %  HERE IS AN INTERESTING WAY TO VISUALIZE THE ABOVE COVERGENCES
  subplot(4,1,1) ;
    image(256 * mvf1') ;
      colormap(gray(256)) ;
          %  THIS COLOR IS CHOSEN FOR EASY PRINTING, FOR VIEWING
          %  ON THE SCREEN, TRY SOME OF THE OTHER VARIATIONS
      title('Good Circulation') ;

  subplot(4,1,2) ;
    image(256 * mvf2') ;
      colormap(gray(256)) ;
      title('Sticky States') ;

  subplot(4,1,3) ;
    image(256 * mvf3') ;
      colormap(gray(256)) ;
      title('Jump Away Fast') ;

  subplot(4,1,4) ;
    image(256 * mvf4') ;
      colormap(gray(256)) ;
      title('Nearly Cyclic') ;


  print s322eg9.ps -dps ;

end ;
