1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
|
## Copyright (C) 2008, 2009, 2010, 2011, 2012, 2016, 2018 Moreno Marzolla
##
## This file is part of the queueing toolbox.
##
## The queueing toolbox is free software: you can redistribute it and/or
## modify it under the terms of the GNU General Public License as
## published by the Free Software Foundation, either version 3 of the
## License, or (at your option) any later version.
##
## The queueing toolbox is distributed in the hope that it will be
## useful, but WITHOUT ANY WARRANTY; without even the implied warranty
## of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with the queueing toolbox. If not, see <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
##
## @deftypefn {Function File} {@var{t} =} ctmcmtta (@var{Q}, @var{p})
##
## @cindex Markov chain, continuous time
## @cindex continuous time Markov chain
## @cindex CTMC
## @cindex mean time to absorption, CTMC
##
## Compute the Mean-Time to Absorption (MTTA) of the CTMC described by
## the infinitesimal generator matrix @var{Q}, starting from initial
## occupancy probabilities @var{p}. If there are no absorbing states, this
## function fails with an error.
##
## @strong{INPUTS}
##
## @table @code
##
## @item @var{Q}(i,j)
## @math{N \times N} infinitesimal generator matrix. @code{@var{Q}(i,j)}
## is the transition rate from state @math{i} to state @math{j}, @math{i
## \neq j}. The matrix @var{Q} must satisfy the condition
## @math{\sum_{j=1}^N Q_{i,j} = 0}
##
## @item @var{p}(i)
## probability that the system is in state @math{i}
## at time 0, for each @math{i=1, @dots{}, N}
##
## @end table
##
## @strong{OUTPUTS}
##
## @table @code
##
## @item @var{t}
## Mean time to absorption of the process represented by matrix @var{Q}.
## If there are no absorbing states, this function fails.
##
## @end table
##
## @strong{REFERENCES}
##
## @itemize
## @item
## G. Bolch, S. Greiner, H. de Meer and
## K. Trivedi, @cite{Queueing Networks and Markov Chains: Modeling and
## Performance Evaluation with Computer Science Applications}, Wiley,
## 1998.
## @end itemize
##
## @seealso{ctmcexps}
##
## @end deftypefn
## Author: Moreno Marzolla <moreno.marzolla(at)unibo.it>
## Web: http://www.moreno.marzolla.name/
function t = ctmcmtta( Q, p )
persistent epsilon = 10*eps;
if ( nargin != 2 )
print_usage();
endif
[N err] = ctmcchkQ(Q);
(N>0) || ...
error(err);
( isvector(p) && length(p) == N && all(p>=0) && abs(sum(p)-1.0)<epsilon ) || ...
error( "p must be a probability vector" );
p = p(:)';
L = ctmcexps(Q,p);
t = sum(L);
endfunction
%!test
%! Q = [0 1 0; 1 0 1; 0 1 0 ]; Q -= diag( sum(Q,2) );
%! fail( "ctmcmtta(Q,[1 0 0])", "no absorbing");
%!test
%! Q = [0 1 0; 1 0 1; 0 0 0; 0 0 0 ];
%! fail( "ctmcmtta(Q,[1 0 0])", "square matrix");
%!test
%! Q = [0 1 0; 1 0 1; 0 0 0 ];
%! fail( "ctmcmtta(Q,[1 0 0])", "infinitesimal");
%!test
%! Q = [ 0 0.1 0 0; ...
%! 0.9 0 0.1 0; ...
%! 0 0.9 0 0.1; ...
%! 0 0 0 0 ];
%! Q -= diag( sum(Q,2) );
%! assert( ctmcmtta( Q,[0 0 0 1] ), 0 ); # state 4 is absorbing
%!test
%! Q = [-1 1; 0 0];
%! assert( ctmcmtta( Q, [0 1] ), 0 ); # state 2 is absorbing
%! assert( ctmcmtta( Q, [1 0] ), 1 ); # the result has been computed by hand
## Compute the MTTA of a pure death process with 4 states
## (state 1 is absorbing). State 4 is the initial state.
%!demo
%! mu = 0.01;
%! death = [ 3 4 5 ] * mu;
%! birth = 0*death;
%! Q = ctmcbd(birth,death);
%! t = ctmcmtta(Q,[0 0 0 1])
%!demo
%! N = 100;
%! birth = death = ones(1,N-1); birth(1) = death(N-1) = 0;
%! Q = diag(birth,1)+diag(death,-1);
%! Q -= diag(sum(Q,2));
%! t = zeros(1,N/2);
%! initial_state = 1:(N/2);
%! for i=initial_state
%! p = zeros(1,N); p(i) = 1;
%! t(i) = ctmcmtta(Q,p);
%! endfor
%! plot(initial_state,t,"+");
%! xlabel("Initial state");
%! ylabel("MTTA");
|
