% 
% Written by:
% -- 
% John L. Weatherwax                2007-09-10
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

close all; drawnow; 
clc; clear; 

format short; 

warning('off'); 
% mc_sample.m
addpath( '~/Software/FullBNT-1.0.4/HMM' ); 
addpath('~/Programming/Matlab/Code/ExternalSources/Engineering/BayesianNetworks/FullBNT/hmm');
% sample_discrete.m
addpath( '~/Software/FullBNT-1.0.4/KPMtools' ); 
% cshistden.m
addpath( '~/Martinez/CSHWM/Code/CSTool' ); 
warning('on'); 

rand('seed',0); randn('seed',0); 

% our transformation matrix in natural form is: 
% -- state i means that we have i normal subgenes from 3
% 
P = [   1,   0,   0,   0; 
      1/5, 3/5, 1/5,   0; 
        0, 1/5, 3/5, 1/5; 
        0,   0,   0,   1  
    ]; 

nTraj = 100000; 
prior    = [ 0, 0, 1, 0 ]; % <- we have two normal genes 
%prior    = [ 0, 0.5, 0.5, 0 ]; % <- equaily in state 1 or 2 

nTimeSteps = 30; 

% generate nTraj Markov chain samples: 
samples = mc_sample( prior, P, nTimeSteps, nTraj ); 

% the number of timesteps we take needs to be large enough so that we end in a
% absorbing state.  We can evaluate empirically if we have selected a large
% enough number with the following code
pNotFinished = sum( ~( (samples(:,end)==4) | (samples(:,end)==1) ) )/nTraj; 
pFinished = 1 - pNotFinished; 
fprintf('an estimate of the probability we arrived at an absorbing state for\n');
fprintf('all the chains in our sample is %10.6f\n', pFinished);


% Part (a):
%
% count the number that end in state "3" (or 4 in ones based indexing)
%
p_3 = sum(samples(:,end)==4)/nTraj; 
fprintf('estimate of the probability we end with all normal subgenes is %10.6f\n', p_3); 
fprintf('the true value is %10.6f diff = %10.6f\n', 2/3, 2/3 - p_3); 

% Part (b):
% 
% count the number of generations before we reach an absorbing state:
%
num_of_trans = zeros(1,nTraj); 
for ii=1:nTraj, 
  inds = find( (samples(ii,:)==1) | (samples(ii,:)==4) );
  if( isempty(inds) )
    num_of_trans(ii) = nTimeSteps+1; % <- this is a lower bound on the number of steps we must take to reach an absorbing state
  else
    num_of_trans(ii) = inds(1);   
  end
end
num_of_trans = num_of_trans-1;       % <- compute the number of transitions from the index 

fprintf('the average of the time to reach an absorbing state is %10.6f\n', mean(num_of_trans));
fprintf('the true value is %10.6f diff = %10.6f\n', 5, 5 - mean(num_of_trans)); 
fprintf('the standard deviation of the time to reach an absorbing state is %10.6f\n', std(num_of_trans));

%figure; hist( num_of_trans ); 
if( 0 ) 
  figure; cshistden( num_of_trans );
  xlabel( 'the number of timesteps until we reach a terminating state' );
  ylabel( 'probability density' ); 
  saveas( gcf, '../../WriteUp/Graphics/Chapter6/chap_6_prob_26_term_state_dist', 'epsc' ); 
end