%
%  macm 202 -- 04 feb 03 -- djm
%
%  w5exp.m:  simulate exponential transition process
%
%     - use state diagram of figure 6.1

clear
mu = 1.0;

%  bins for histogram & pseudorandom seed

dt = 1/8;  mudt = mu*dt;
bins = (0:dt/mu:5/mu);
%rand('state',7654321);

%  simulation of Nruns at dt

Nruns = 10000;
wait_times = zeros(Nruns,1);

%  simulation loop of figure 6.1

for k=1:Nruns;
    wait = 0;
    while (rand(1,1)>mudt)
        wait = wait + dt;
    end
    wait_times(k,1) = wait;
end    

%  plot histogram of transition probabilities

figure(1); clf
count = hist(wait_times,bins);
bar(bins,count/Nruns,'hist')
hold on
plot(bins,exp(-mu*bins)*dt,'r')
plot(bins,(1-mudt).^(bins/dt)*dt,'g')
title('\bf simulation of exponential transition probability')
xlabel('\bf time, t')
ylabel('\bf p(transition during t->t+\Delta t)')
legend('\bf exponential limit','\bf exact discrete solution')
axis([0 5 0 1.1*dt])

%  direct generation of Nruns exponential transition times!

rand_wait = -log(rand(Nruns,1))/mu;

figure(2); clf
count = histc(rand_wait,bins);
bar(bins,count/Nruns,'histc')
hold on
plot(bins,exp(-mu*bins)*dt,'r')
title('\bf direct generation of exponential transition probability')
xlabel('\bf time, t')
ylabel('\bf p(transition during t->t+\Delta t)')
legend('\bf exponential limit')
axis([0 5 0 1.1*dt])