% This code solves problem 13 of chapter 4. The correct answer is
% approximately 2.56.
%author: David Anderson
%Usage:
% [ExpectedTimeOfDeath,TwoStandardDeviations] =
% BranchingExtinction(10000);
%
function [ExpectedTimeOfDeath,TwoStandardDeviations] = BranchingExtinction(n)
%Maximum number of steps to take in the branching process.
maxi = 10^5;
%collect the data.
data = zeros([n,1]);
%Perform n experiments for the purposes of Monte Carlo.
for k = 1:n
%Reset population to zero.
X = zeros([maxi,1]);
%Initial population to one.
X(1) = 1;
%I will use this to store the extinction times.
Value = 0;
%Generate a trajectory of the branching process.
for i =1:maxi-1
%If the process went to zero, return the time (later I will
%subtract off one since we started at time one)
if X(i) == 0
Value = i;
break
end
%I need a vector holding the number of offspring of the total
%number of people in the population.
Y = zeros([X(i),1]);
%we have to find out how many offspring there were for each person
% in the population.
for j = 1:X(i)
r = rand;
if r < 1/2
Y(j) = 0;
elseif r < 3/4
Y(j) = 1;
else
Y(j) = 2;
end
end
%Update the population size.
X(i+1)= sum(Y);
end
%If I want to plot realizations, allow the code below (but make sure to
%not have a large value of n.
% Tn = zeros([Value,1]);
% for i = 1:Value-1
% Tn(i+1) = Tn(i) + 1;
% end
%
% size(X(1:Value))
% size(Tn)
%
% figure
% plot(Tn,X(1:Value))
%subtract off one to account for the fact that we started at time 1
%instead of time zero.
Value = Value - 1;
%update data vector.
data(k) = Value;
end
%Find mean of data
ExpectedTimeOfDeath=mean(data)
%find two standard deviations of the data.
TwoStandardDeviations = 1.96*sqrt(var(data)/n)
end %The program