% 
% 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.
% 
%-----

clc; 
close all; drawnow; 
clear all; 

% some constants of our birth death process: 
% 
lambda0 = 1; 
mu0     = 2; 
M       = 4; 
s       = 2; 

% the time to simulate to: 
t_stop = 1.0; 

% the current time (and an array of event times): 
t = 0; 
t_arr = [t]; 
% the current state (and a array of state transitions):
n = 0; 
n_arr = [n]; 
while( 1 )
  if( t >= t_stop ) break; end 
  % get the time till the next transformation: 
  soj_rate = chap_9_prob_3_lambda(n,M,lambda0) + chap_9_prob_3_mu(n,s,mu0); 
  T = exprnd( 1/soj_rate, 1, 1 ); 
  t_arr(end+1) = T; 
  % determine the type of transformation it was (and update the state): 
  p_nP1 = chap_9_prob_3_lambda(n,M,lambda0)/soj_rate; 
  if( rand(1) < p_nP1 )
    n=n+1; 
  else
    n=n-1; 
  end
  n_arr(end+1) = n; 
  t = t+T; 
end

t_cs = cumsum(t_arr); 

v = [ t_cs(:), t_cs(:) ].'; 
v = v(:); v(1)=[]; 
y = [ n_arr(:), n_arr(:) ].';
y = y(:); y(end)=[]; 

figure; plot( v, y, '-x', 'markersize', 5 ); 
xlabel( 't (time)' ); ylabel( 'X(t)' ); axis( [min(t_cs),t_stop,0,max(n_arr)+1] ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter9/chap_9_prob_3_xt', 'epsc' );