Simulation of Markov chains

Question

I have the following Markov chain:

This chain shows the states of the Spaceship, which is in the asteroid belt: S1 - is serviceable, S2 - is broken. 0.12 - the probability of destroying the Spaceship by a collision with an asteroid. 0.88 - the probability of that a collision will not be critical. Need to find the probability of a serviceable condition of the ship after the third collision.

Analytical solution showed the response - 0.681. But it is necessary to solve this problem by simulation method using any modeling tool (MATLAB Simulink, AnyLogic, Scilab, etc.).

Do you know what components should be used to simulate this process in Simulink or any other simulation environment? Any examples or links.

Answer 1

For a more complex system you'll want to use Stateflow or SimEvents, but for this simple example all you need is a single Unit Delay block (output = 0 => S1, output = 1 => S2), with a Switch block, a Random block, and some comparison blocks to construct the logic determining the next value of the state.

Presumably you must execute the simulation a (very) large number of times and average the results to get a statistically significant output. You'll need to change the "seed" of the random generator each time you run the simulation. This can be done by setting the seed to be "now" (or something similar to that).

Alternatively you could quite easily vectorize the model so that you only need to execute it once.

Answer 2

First, we know the three step probability transition matrix contains the answer (0.6815).

% MATLAB R2019a
P = [0.88 0.12;
    0 1];
P3 = P*P*P
P(1,1)                 % 0.6815

Approach 1: Requires Econometrics Toolbox
This approach uses the dtmc() and simulate() functions.

First, create the Discrete Time Markov Chain (DTMC) with the probability transition matrix, P, and using dtmc().

mc = dtmc(P);               % Create the DTMC
numSteps = 3;               % Number of collisions

You can get one sample path easily using simulate(). Pay attention to how you specify the initial conditions.

% One Sample Path
rng(8675309)                                     % for reproducibility
X = simulate(mc,numSteps,'X0',[1 0])

% Multiple Sample Paths
numSamplePaths = 3;
X = simulate(mc,numSteps,'X0',[numSamplePaths 0])  % returns a 4 x 3 matrix

The first row is the X0 row for the starting state (initial condition) of the DTMC. The second row is the state after 1 transition (X1). Thus, the fourth row is the state after 3 transitions (collisions).

% 50000 Sample Paths
rng(8675309)                                     % for reproducibility
k = 50000;
X = simulate(mc,numSteps,'X0',[k 0]);            % returns a 4 x 50000 matrix

prob_survive_3collisions = sum(X(end,:)==1)/k    % 0.6800

We can bootstrap a 95% Confidence Interval on the mean probability to survive 3 collisions to get 0.6814 ± 0.00069221, or rather, [0.6807 0.6821], which contains the result.

numTrials = 40;
ProbSurvive_3collisions = zeros(numTrials,1);
for trial = 1:numTrials
    Xtrial = simulate(mc,numSteps,'X0',[k 0]);
    ProbSurvive_3collisions(trial) = sum(Xtrial(end,:)==1)/k;
end

% Mean +/- Halfwidth
alpha = 0.05;
mean_prob_survive_3collisions = mean(ProbSurvive_3collisions)
hw = tinv(1-(0.5*alpha), numTrials-1)*(std(ProbSurvive_3collisions)/sqrt(numTrials))
ci95 = [mean_prob_survive_3collisions-hw mean_prob_survive_3collisions+hw]   

---

maxNumCollisions = 10;
numSamplePaths = 50000;
ProbSurvive = zeros(maxNumCollisions,1);
for numCollisions = 1:maxNumCollisions
    Xc = simulate(mc,numCollisions,'X0',[numSamplePaths 0]);
    ProbSurvive(numCollisions) = sum(Xc(end,:)==1)/numSamplePaths;
end

Answer 3

If you want to simulate this, it is fairly easy in matlab:

servicable = 1;
t = 0;
while servicable =1
   t = t+1;
  servicable = rand()<=0.88
end

Now t represents the amount of steps before the ship is broken.

Wrap this in a for loop and you can do as many simulations as you like.

---

Note that this can actually give you the distribution, if you want to know it after 3 times, simply add && t<3 to the while condition.