disp(['Running MATLAB script s322hw7.m']) ;
%
%    SOLUTION TO HOMEWORK 6 FOR STATISTICS 322
%    POSTED 10/24/96
%    EXPLORES MONTE CARLO MARKOV CHAIN

format compact ;

ipart = 2 ;      %  1,2

if ipart == 1 ;    %  THEN MAKE PLOTS OF SINGLE PATHS

  rand('seed',87309385) ;

  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 ;
          %  THIS BLOCK COPIED DIRECTLY FROM s322eg9.m



    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]


    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


    x2 = 1 ;
          %  ARBITRARY STARTING VALUE

    mx = [] ;
    nstep = 100 ;
    for istep = 1:nstep ;
      vp = f21(:,x2) ;
          %  THE X2-TH COLUMN GIVES THE VECTOR OF PROBS, GIVEN X2
      unif = rand(1) ;
      if unif <= vp(1) ;
        x1 = 1 ;
      elseif unif <= vp(1) + vp(2) ;
        x1 = 2;
      else ;
        x1 = 3;
      end ;

      vp = f12(:,x1) ;
          %  THE X1-TH COLUMN GIVES THE VECTOR OF PROBS, GIVEN X1
      unif = rand(1) ;
      if unif <= vp(1) ;
        x2 = 1 ;
      elseif unif <= vp(1) + vp(2) ;
        x2 = 2;
      else ;
        x2 = 3;
      end ;

      mx = [mx; x1 x2] ;
    end ;

    subplot(4,1,ijoint) ;
      plot(1:nstep,mx(:,1),'-',1:nstep,mx(:,2),':') ;
        title(fstr) ;
        axis([0, 100, 0, 4]) ;

  end ;

  print s322hw7.ps -dps ;


elseif ipart == 2 ;    %  THEN DO CHECK OF MONTE CARLO VARIABILITY

  rand('seed',34586351) ;

  %  WRITE ANSWERS INTO A FILE
  fid = fopen('s322hw7.out','wt') ;
    titlstr1 = 'Output from the MATLAB Script s322hw7.m,   ' ;
    titlstr1 = [titlstr1,date] ;
  cntbytes = fprintf(fid,'%1s\n',titlstr1) ;
    titlstr2 = '      Simulations for Monte Carlo Markov Chains, ' ;
  cntbytes = fprintf(fid,'%1s\n\n\n',titlstr2) ;
    titlstr3 = '   joint dist       sd(xbar)  bias(xbar)     ' ;
    titlstr3 = [titlstr3 'sc(cov)   bias(cov)'] ;
  cntbytes = fprintf(fid,'%1s\n',titlstr3) ;

  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 ;
          %  THIS BLOCK COPIED DIRECTLY FROM s322eg9.m


    disp(['For ' fstr]) ;

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

    ex1 = f1' * (1:3)' ;
    ex2 = f2 * (1:3)' ;
    ex1x2 = (1:3) * f * (1:3)' ;
    tcov = ex1x2 - ex1 * ex2 ;

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


    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


    nsim = 200 ;
    vx2 = ones(1,nsim) ;
          %  ARBITRARY STARTING VALUE

    mx1 = [] ;
    mx2 = [] ;
    nburnin = 40 ;
    nstep = 400 ;
    for istep = 1:(nburnin+nstep) ;
      vp = f21(:,vx2) ;
          %  THE X2-TH COLUMN GIVES THE VECTOR OF PROBS, GIVEN X2
      unif = rand(1,nsim) ;
      vx1 = 2 * ones(1,nsim) ;
          %  START WITH VECTOR OF 2'S
        flag = find(unif <= vp(1,:)) ;
      vx1(1,flag) = ones(1,length(flag)) ;
          %  REPLACE SOME WITH 1'S
        flag = find(unif > (1 - vp(3,:))) ;
      vx1(1,flag) = 3 * ones(1,length(flag)) ;
          %  REPLACE OTHERS WITH 3'S

      vp = f12(:,vx1) ;
          %  THE X1-TH COLUMN GIVES THE VECTOR OF PROBS, GIVEN X1
      unif = rand(1,nsim) ;
      vx2 = 2 * ones(1,nsim) ;
          %  START WITH VECTOR OF 2'S
        flag = find(unif <= vp(1,:)) ;
      vx2(1,flag) = ones(1,length(flag)) ;
          %  REPLACE SOME WITH 1'S
        flag = find(unif > (1 - vp(3,:))) ;
      vx2(1,flag) = 3 * ones(1,length(flag)) ;
          %  REPLACE OTHERS WITH 3'S

      mx1 = [mx1; vx1] ;
      mx2 = [mx2; vx2] ;
          %  CREATE BIG MATRICES OF DATA
    end ;

    mx1 = mx1((nburnin+1):(nburnin+nstep),:) ;
    mx2 = mx2((nburnin+1):(nburnin+nstep),:) ;
          %  KEEP ONLY PART AFTER BURN IN

    vx1bar = mean(mx1) ; 
    x1bar = mean(vx1bar) ;
    sdxb = std(vx1bar) 
    biasxb = x1bar - ex1
          %  WRITE THESE TWO TO SCREEN

    vx2bar = mean(mx2) ;
    vx1x2bar = mean(mx1 .* mx2) ;
    vscov = vx1x2bar - vx1bar .* vx2bar ;
    sdsc = std(vscov)
    biassc = mean(vscov) - tcov 
          %  WRITE THESE TWO TO SCREEN
 

    %  WRITE TO FILE
      forstr = [fstr '    %6.3f     (%6.3f)      %6.3f    (%6.3f)\n'] ;
      vout = [sdxb biasxb sdsc biassc] ;
    cntbytes = fprintf(fid,forstr,vout) ;

  end ;

  fclose('all') ;

end ; 

