| << Chapter < Page | Chapter >> Page > |
branchp.m Calculates the transition matrix for a simple branching process with a specified maximum population. Input consists of the maximum population value M and the coefficient matrix for the generating function for the individual propagation random variables Z i . The latter matrix must include zero coefficients for missing powers.
% BRANCHP file branchp.m Transition P for simple branching process
% Version of 7/25/95% Calculates transition matrix for a simple branching
% process with specified maximum population.disp('Do not forget zero probabilities for missing values of Z')
PZ = input('Enter PROBABILITIES for individuals ');M = input('Enter maximum allowable population ');
mz = length(PZ) - 1;EZ = dot(0:mz,PZ);
disp(['The average individual propagation is ',num2str(EZ),])
P = zeros(M+1,M+1);Z = zeros(M,M*mz+1);
k = 0:M*mz;a = min(M,k);
z = 1;P(1,1) = 1;
for i = 1:M % Operation similar to gendz = conv(PZ,z);
Z(i,1:i*mz+1) = z;[t,p] = csort(a,Z(i,:));P(i+1,:) = p;
enddisp('The transition matrix is P')
disp('To study the evolution of the process, call for branchdbn') chainset.m Sets up for simulation of Markov chains. Inputs are the transition matrix P the set of states, and an optional set of target states. The chain generating procedures listed below assume this procedure has been run.
% CHAINSET file chainset.m Setup for simulating Markov chains
% Version of 1/2/96 Revise 7/31/97 for version 4.2 and 5.1P = input('Enter the transition matrix ');
ms = length(P(1,:));MS = 1:ms;
states = input('Enter the states if not 1:ms ');if isempty(states)
states = MS;end
disp('States are')disp([MS;states]')PI = input('Enter the long-run probabilities ');
F = [zeros(1,ms); cumsum(P')]';
A = F(:,MS);B = F(:,MS+1);
e = input('Enter the set of target states ');ne = length(e);
E = zeros(1,ne);for i = 1:ne
E(i) = MS(e(i)==states);end
disp(' ')disp('Call for for appropriate chain generating procedure') mchain.m Assumes chainset has been run. Generates trajectory of specified length, with specified initial state.
% MCHAIN file mchain.m Simulation of Markov chains
% Version of 1/2/96 Revised 7/31/97 for version 4.2 and 5.1% Assumes the procedure chainset has been run
n = input('Enter the number n of stages ');st = input('Enter the initial state ');
if ~isempty(st)s = MS(st==states);
elses = 1;
endT = zeros(1,n); % Trajectory in state numbers
U = rand(1,n);for i = 1:n
T(i) = s;s = ((A(s,:)<U(i))&(U(i)<= B(s,:)))*MS';
endN = 0:n-1;
tr = [N;states(T)]';
n10 = min(n,11);TR = tr(1:n10,:);
f = ones(1,n)/n;[sn,p] = csort(T,f);if isempty(PI)
disp(' State Frac')disp([states; p]')else
disp(' State Frac PI')disp([states; p; PI]')end
disp('To view the first part of the trajectory of states, call for TR') arrival.m Assumes chainset has been run. Calculates repeatedly the arrival time to a prescribed set of states.
% ARRIVAL file arrival.m Arrival time to a set of states
% Version of 1/2/96 Revised 7/31/97 for version 4.2 and 5.1% Calculates repeatedly the arrival
% time to a prescribed set of states.% Assumes the procedure chainset has been run.
r = input('Enter the number of repetitions ');disp('The target state set is:')
disp(e)st = input('Enter the initial state ');
if ~isempty(st)s1 = MS(st==states); % Initial state number
elses1 = 1;
endclear T % Trajectory in state numbers (reset)
S = zeros(1,r); % Arrival time for each rep (reset)TS = zeros(1,r); % Terminal state number for each rep (reset)
for k = 1:rR = zeros(1,ms); % Indicator for target state numbers
R(E) = ones(1,ne); % reset for target state numberss = s1;
T(1) = s;i = 1;
while R(s) ~= 1 % While s is not a target state numberu = rand(1,1);
s = ((A(s,:)<u)&(u<= B(s,:)))*MS';
i = i+1;T(i) = s;
endS(k) = i-1; % i is the number of stages; i-1 is time
TS(k) = T(i);end
[ts,ft]= csort(TS,ones(1,r)); % ts = terminal state numbers ft = frequencies
fts = ft/r; % Relative frequency of each ts[a,at] = csort(TS,S); % at = arrival time for each tsw = at./ft; % Average arrival time for each ts
RES = [states(ts); fts; w]';
disp(' ')if r == 1
disp(['The arrival time is ',int2str(i-1),])
disp(['The state reached is ',num2str(states(ts)),])
N = 0:i-1;TR = [N;states(T)]';disp('To view the trajectory of states, call for TR')
elsedisp(['The result of ',int2str(r),' repetitions is:'])disp('Term state Rel Freq Av time')
disp(RES)disp(' ')
[t,f]= csort(S,ones(1,r)); % t = arrival times f = frequencies
p = f/r; % Relative frequency of each tdbn = [t; p]';AV = dot(t,p);
SD = sqrt(dot(t.^2,p) - AV^2);MN = min(t);
MX = max(t);disp(['The average arrival time is ',num2str(AV),])disp(['The standard deviation is ',num2str(SD),])disp(['The minimum arrival time is ',int2str(MN),])disp(['The maximum arrival time is ',int2str(MX),])disp('To view the distribution of arrival times, call for dbn')
disp('To plot the arrival time distribution, call for plotdbn')end recurrence.m Assumes chainset has been run. Calculates repeatedly the recurrence time to a prescribed set of states, if initial state is in the set; otherwise calculatesthe arrival time.
% RECURRENCE file recurrence.m Recurrence time to a set of states
% Version of 1/2/96 Revised 7/31/97 for version 4.2 and 5.1% Calculates repeatedly the recurrence time
% to a prescribed set of states, if initial% state is in the set; otherwise arrival time.
% Assumes the procedure chainset has been run.r = input('Enter the number of repititions ');
disp('The target state set is:')disp(e)
st = input('Enter the initial state ');if ~isempty(st)
s1 = MS(st==states); % Initial state numberelse
s1 = 1;end
clear T % Trajectory in state numbers (reset)S = zeros(1,r); % Recurrence time for each rep (reset)
TS = zeros(1,r); % Terminal state number for each rep (reset)for k = 1:r
R = zeros(1,ms); % Indicator for target state numbersR(E) = ones(1,ne); % reset for target state numbers
s = s1;T(1) = s;
i = 1;if R(s) == 1
u = rand(1,1);s = ((A(s,:)<u)&(u<= B(s,:)))*MS';
i = i+1;T(i) = s;
endwhile R(s) ~= 1 % While s is not a target state number
u = rand(1,1);s = ((A(s,:)<u)&(u<= B(s,:)))*MS';
i = i+1;T(i) = s;
endS(k) = i-1; % i is the number of stages; i-1 is time
TS(k) = T(i);end
[ts,ft]= csort(TS,ones(1,r)); % ts = terminal state numbers ft = frequencies
fts = ft/r; % Relative frequency of each ts[a,tt] = csort(TS,S); % tt = total time for each tsw = tt./ft; % Average time for each ts
RES = [states(ts); fts; w]';
disp(' ')if r == 1
disp(['The recurrence time is ',int2str(i-1),])
disp(['The state reached is ',num2str(states(ts)),])
N = 0:i-1;TR = [N;states(T)]';disp('To view the trajectory of state numbers, call for TR')
elsedisp(['The result of ',int2str(r),' repetitions is:'])disp('Term state Rel Freq Av time')
disp(RES)disp(' ')
[t,f]= csort(S,ones(1,r)); % t = recurrence times f = frequencies
p = f/r; % Relative frequency of each tdbn = [t; p]';AV = dot(t,p);
SD = sqrt(dot(t.^2,p) - AV^2);MN = min(t);
MX = max(t);
disp(['The average recurrence time is ',num2str(AV),])
disp(['The standard deviation is ',num2str(SD),])
disp(['The minimum recurrence time is ',int2str(MN),])
disp(['The maximum recurrence time is ',int2str(MX),])
disp('To view the distribution of recurrence times, call for dbn')disp('To plot the recurrence time distribution, call for plotdbn')
end kvis.m Assumes chainset has been run. Calculates repeatedly the time to complete visits to a specified k of the states in a prescribed set.
% KVIS file kvis.m Time to complete k visits to a set of states
% Version of 1/2/96 Revised 7/31/97 for version 4.2 and 5.1% Calculates repeatedly the time to complete
% visits to k of the states in a prescribed set.% Default is visit to all the target states.
% Assumes the procedure chainset has been run.r = input('Enter the number of repetitions ');
disp('The target state set is:')disp(e)
ks = input('Enter the number of target states to visit ');if isempty(ks)
ks = ne;end
if ks>ne
ks = ne;end
st = input('Enter the initial state ');if ~isempty(st)
s1 = MS(st==states); % Initial state numberelse
s1 = 1;end
disp(' ')clear T % Trajectory in state numbers (reset)
R0 = zeros(1,ms); % Indicator for target state numbersR0(E) = ones(1,ne); % reset
S = zeros(1,r); % Terminal transitions for each rep (reset)for k = 1:r
R = R0;s = s1;
if R(s) == 1R(s) = 0;
endi = 1;
T(1) = s;while sum(R)>ne - ks
u = rand(1,1);s = ((A(s,:)<u)&(u<= B(s,:)))*MS';
if R(s) == 1R(s) = 0;
endi = i+1;
T(i) = s;end
S(k) = i-1;end
if r == 1disp(['The time for completion is ',int2str(i-1),])N = 0:i-1;
TR = [N;states(T)]';
disp('To view the trajectory of states, call for TR')else
[t,f]= csort(S,ones(1,r));
p = f/r;D = [t;f]';AV = dot(t,p);
SD = sqrt(dot(t.^2,p) - AV^2);MN = min(t);
MX = max(t);disp(['The average completion time is ',num2str(AV),])disp(['The standard deviation is ',num2str(SD),])disp(['The minimum completion time is ',int2str(MN),])disp(['The maximum completion time is ',int2str(MX),])disp(' ')
disp('To view a detailed count, call for D.')disp('The first column shows the various completion times;')
disp('the second column shows the numbers of trials yielding those times')end plotdbn Used after m-procedures arrival or recurrence to plot arrival or recurrence time distribution.
% PLOTDBN file plotdbn.m
% Version of 1/23/98% Plot arrival or recurrence time dbn
% Use after procedures arrival or recurrence% to plot arrival or recurrence time distribution
plot(t,p,'-',t,p,'+')grid
title('Time Distribution')xlabel('Time in number of transitions')
ylabel('Relative frequency') 
Notification Switch
Would you like to follow the 'Applied probability' conversation and receive update notifications?
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|