In our hack session I coded a small CTMC simulator in Matlab to implement the 6-state example from Hao's group. The goal of this is to have a starting case example to develop importance sampling approaches. We ended with some confusion about how to check "mean time to fail" in PRISM. It turns out to be fairly simple using reward syntax. In this post I have the Matlab simulation code and also the PRISM version, and confirm that they get similar results (I say "similar" because I haven't done a rigorous statistical comparison, just ta few cases).
Here is the Matlab simulation:
% Example Markov Process
% has 6 states
states = 0:6;
labels = {
'init'
''
'busy'
''
''
''
'crash'
};
% Rate matrix:
% 0 1 2 3 4 5 6
R = [0 8 0 2 0 0 0 % 0
0 0 8 0 1 0 0 % 1
0 2 0 0 0 1 0 % 2
0 0 0 0 0 0 9 % 3
0 0 0 0 0 0 9 % 4
0 0 0 0 0 0 9 % 5
0 0 0 0 0 0 0 % 6
];
% initial state
s0 = 0;
% N = number of paths to simulate
% Each path is run until we reach a "target"
N = 10000;
simulation_history = {};
mean_path_length = 0;
% Simulate the paths
for path=1:N
target = 0;
s = s0;
t = 0;
trace_index = 0;
path_trace = [t s]';
% This simulation terminates paths when they reach
% the "crash" state, so it can be used to compute
% (estimate) mean-time-to-crash
while (~target)
% In a CTMC, the time-to-transition is
% exponentially distributed
possible_transitions = find(R(s+1,:)>0);
rates = R(s+1,possible_transitions);
% Generate some exponential random numbers:
r = exprnd(1./rates);
% Pick the transition that happens soonest:
actual_transition = possible_transitions(find(r==min(r)));
% Advance the simulation time and state:
t = t + min(r);
s = actual_transition-1;
% Add this transition to the simulation trace:
trace_index = trace_index + 1;
path_trace(:,trace_index) = [t s]';
% Check termination condition:
if (s == 6)
target = 1;
end
%fprintf("Time %f\t State %d\n",t,s);
%pause
end
% Add path trace to ensemble:
simulation_history{path} = path_trace;
% Add to mean_path_length:
mean_path_length = mean_path_length + t/N;
end
fprintf('Mean path duration: %f\n',mean_path_length)
% What we really want is the probability of violating a
% given property. In this example, that property is "true U<p crash"
% P<0.22 [ !"busy" U "crash" ] "(not busy) until crash", i.e.
% crashing without ever being busy, so we don't care about any paths
% that traverse state 2.
% For importance sampling:
%
% 1. Suppress part of the graph that contradicts the property
% 2. Enhance part of the graph indicated by counter-examples
The corresponding PRISM model is shown below. This is identical to the ``example.sm`` model except a reward is introduced where every non-crash state gets a reward of "1", so the total reward should be the total path length.
ctmc
module test
s : [0..6] init 0;
[] s=0 -> 8.0 : (s'=1) + 2.0 : (s'=3);
[] s=1 -> 8.0 : (s'=2) + 1.0 : (s'=4);
[] s=2 -> 2.0 : (s'=1) + 1.0 : (s'=5);
[] s=3 -> 9.0 : (s'=6);
[] s=4 -> 9.0 : (s'=6);
[] s=5 -> 9.0 : (s'=6);
endmodule
label "ready" = s=1 ;
label "busy" = s=2 ;
label "crash" = s=6 ;
rewards "total_time"
s<6 : 1;
endrewards
To verify the mean path length:
R=? [ C ]The Matlab and PRISM implementations agree with a mean path length of 1.01. When changing the transition rates, e.g. changing the rate at (1,4) from "1" to "4" gives a mean path length of 0.65 in both platforms.
#
examples #
Matlab #
PRISM #
CTMC #
simulation #
MTTF #
HackSession
