%
%  math495/stat490 -- 16 nov 03 -- djm
%
%  w12ehr.m:  

clear

%  test parameters
Nexpt = 10000;  Nitems = 4;

%  initialize random number generator
phone = 2914814;
rand('state',phone);  %  initialize with my office phone number

%  data
data = zeros(Nexpt,1);

%  run ehrenfest experiments
for kk=1:Nexpt
    items = 0;  moves = 0;
    while (items~=Nitems)
        items = items - 1 + 2*(rand(1,1)>items/Nitems);
        moves = moves + 1;
    end
    data(kk) = moves;
end

figure(1)
count = histc(data,2:2:102);
subplot(2,1,1)
plot(2:2:100,count(1:end-1),'*')
axis([0 100 0 Nexpt/10])
title('\bf DPF of first exit times')
xlabel('\bf exit times')

subplot(2,1,2)
plot(4:2:100,log(count(2:end-1)),'*')
axis([0 100 0 10])
[const] = polyfit(4:2:100,log(count(2:end-1))',1);
slope = const(1)
title('\bf semi-log plot:  log|DPF(n)| vs n')
ylabel('\bf log(DPF)')
xlabel('\bf exit times')
text(40,8,['\bf slope = ' num2str(slope)])

%  4-item expectations
if (Nitems==4)
    Q = [0 1 0 0 ; 1/4 0 3/4 0 ; 0 1/2 0 1/2 ; 0 0 3/4 0];
    b = ones(4,1);

    E = (eye(4,4)-Q)\b;
    V = (eye(4,4)-Q)\(2*E-b) - E.*E;
end

exit_mean = [mean(data) E(1)]
exit_var  = [var(data) V(1)]