1## Copyright (C) 2011, 2012, 2016, 2018 Moreno Marzolla 2## 3## This file is part of the queueing toolbox. 4## 5## The queueing toolbox is free software: you can redistribute it and/or 6## modify it under the terms of the GNU General Public License as 7## published by the Free Software Foundation, either version 3 of the 8## License, or (at your option) any later version. 9## 10## The queueing toolbox is distributed in the hope that it will be 11## useful, but WITHOUT ANY WARRANTY; without even the implied warranty 12## of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 13## General Public License for more details. 14## 15## You should have received a copy of the GNU General Public License 16## along with the queueing toolbox. If not, see <http://www.gnu.org/licenses/>. 17 18## -*- texinfo -*- 19## 20## @deftypefn {Function File} {@var{M} =} ctmcfpt (@var{Q}) 21## @deftypefnx {Function File} {@var{m} =} ctmcfpt (@var{Q}, @var{i}, @var{j}) 22## 23## @cindex first passage times, CTMC 24## @cindex CTMC 25## @cindex continuous time Markov chain 26## @cindex Markov chain, continuous time 27## 28## Compute mean first passage times for an irreducible continuous-time 29## Markov chain. 30## 31## @strong{INPUTS} 32## 33## @table @code 34## 35## @item @var{Q}(i,j) 36## Infinitesimal generator matrix. @var{Q} is a @math{N \times N} 37## square matrix where @code{@var{Q}(i,j)} is the transition rate from 38## state @math{i} to state @math{j}, for @math{1 @leq{} i, j @leq{} N}, 39## @math{i \neq j}. Transition rates must be nonnegative, and 40## @math{\sum_{j=1}^N Q_{i,j} = 0} 41## 42## @item @var{i} 43## Initial state. 44## 45## @item @var{j} 46## Destination state. 47## 48## @end table 49## 50## @strong{OUTPUTS} 51## 52## @table @code 53## 54## @item @var{M}(i,j) 55## average time before state 56## @var{j} is visited for the first time, starting from state @var{i}. 57## We let @code{@var{M}(i,i) = 0}. 58## 59## @item m 60## @var{m} is the average time before state @var{j} is visited for the first 61## time, starting from state @var{i}. 62## 63## @end table 64## 65## @seealso{ctmcmtta} 66## 67## @end deftypefn 68 69## Author: Moreno Marzolla <moreno.marzolla(at)unibo.it> 70## Web: http://www.moreno.marzolla.name/ 71 72function result = ctmcfpt( Q, i, j ) 73 74 persistent epsilon = 10*eps; 75 76 if ( nargin != 1 && nargin != 3 ) 77 print_usage(); 78 endif 79 80 [N err] = ctmcchkQ(Q); 81 82 (N>0) || ... 83 error(err); 84 85 if ( nargin == 1 ) 86 result = zeros(N,N); 87 for j=1:N 88 QQ = Q; 89 QQ(j,:) = 0; # make state j absorbing 90 for i=1:N 91 p0 = zeros(1,N); p0(i) = 1; 92 result(i,j) = ctmcmtta(QQ,p0); 93 endfor 94 endfor 95 else 96 (isscalar(i) && i>=1 && j<=N) || error("i must be an integer in the range 1..%d", N); 97 (isvector(j) && all(j>=1) && all(j<=N)) || error("j must be an integer or vector with elements in 1..%d", N); 98 j = j(:)'; # make j a row vector 99 Q(j,:) = 0; # make state(s) j absorbing 100 p0 = zeros(1,N); p0(i) = 1; 101 result = ctmcmtta(Q,p0); 102 endif 103endfunction 104%!demo 105%! Q = [ -1.0 0.9 0.1; ... 106%! 0.1 -1.0 0.9; ... 107%! 0.9 0.1 -1.0 ]; 108%! M = ctmcfpt(Q) 109%! m = ctmcfpt(Q,1,3) 110 111%!test 112%! N = 10; 113%! Q = reshape(1:N^2,N,N); 114%! Q(1:N+1:end) = 0; 115%! Q -= diag(sum(Q,2)); 116%! M = ctmcfpt(Q); 117%! assert( all(diag(M) < 10*eps) ); 118 119