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} =} dtmcfpt (@var{P}) 21## 22## @cindex first passage times 23## @cindex mean recurrence times 24## @cindex discrete time Markov chain 25## @cindex Markov chain, discrete time 26## @cindex DTMC 27## 28## Compute mean first passage times and mean recurrence times 29## for an irreducible discrete-time Markov chain over the state space 30## @math{@{1, @dots{}, N@}}. 31## 32## @strong{INPUTS} 33## 34## @table @code 35## 36## @item @var{P}(i,j) 37## transition probability from state @math{i} to state @math{j}. 38## @var{P} must be an irreducible stochastic matrix, which means that 39## the sum of each row must be 1 (@math{\sum_{j=1}^N P_{i j} = 1}), 40## and the rank of @var{P} must be @math{N}. 41## 42## @end table 43## 44## @strong{OUTPUTS} 45## 46## @table @code 47## 48## @item @var{M}(i,j) 49## For all @math{1 @leq{} i, j @leq{} N}, @math{i \neq j}, @code{@var{M}(i,j)} is 50## the average number of transitions before state @var{j} is entered 51## for the first time, starting from state @var{i}. 52## @code{@var{M}(i,i)} is the @emph{mean recurrence time} of state 53## @math{i}, and represents the average time needed to return to state 54## @var{i}. 55## 56## @end table 57## 58## @strong{REFERENCES} 59## 60## @itemize 61## @item Grinstead, Charles M.; Snell, J. Laurie (July 62## 1997). @cite{Introduction to Probability}, Ch. 11: Markov 63## Chains. American Mathematical Society. ISBN 978-0821807491. 64## @end itemize 65## 66## @seealso{ctmcfpt} 67## 68## @end deftypefn 69 70## Author: Moreno Marzolla <moreno.marzolla(at)unibo.it> 71## Web: http://www.moreno.marzolla.name/ 72 73function result = dtmcfpt( P ) 74 75 if ( nargin != 1 ) 76 print_usage(); 77 endif 78 79 [N err] = dtmcchkP(P); 80 81 ( N>0 ) || ... 82 error(err); 83 84 if ( any(diag(P) == 1) ) 85 error("Cannot compute first passage times for absorbing chains"); 86 endif 87 88 w = dtmc(P); # steady state probability vector 89 W = repmat(w,N,1); 90 ## Z = (I - P + W)^-1 where W is the matrix where each row is the 91 ## steady-state probability vector for P 92 Z = inv(eye(N)-P+W); 93 ## m_ij = (z_jj - z_ij) / w_j 94 result = (repmat(diag(Z)',N,1) - Z) ./ repmat(w,N,1) + diag(1./w); 95 96endfunction 97%!test 98%! P = [1 1 1; 1 1 1]; 99%! fail( "dtmcfpt(P)" ); 100 101%!test 102%! P = dtmcbd([1 1 1], [0 0 0] ); 103%! fail( "dtmcfpt(P)", "absorbing" ); 104 105%!test 106%! P = [ 0.0 0.9 0.1; ... 107%! 0.1 0.0 0.9; ... 108%! 0.9 0.1 0.0 ]; 109%! p = dtmc(P); 110%! M = dtmcfpt(P); 111%! assert( diag(M)', 1./p, 1e-8 ); 112 113## Example on p. 461 of [GrSn97] 114%!test 115%! P = [ 0 1 0 0 0; ... 116%! .25 .0 .75 0 0; ... 117%! 0 .5 0 .5 0; ... 118%! 0 0 .75 0 .25; ... 119%! 0 0 0 1 0 ]; 120%! M = dtmcfpt(P); 121%! assert( M, [16 1 2.6667 6.3333 21.3333; ... 122%! 15 4 1.6667 5.3333 20.3333; ... 123%! 18.6667 3.6667 2.6667 3.6667 18.6667; ... 124%! 20.3333 5.3333 1.6667 4 15; ... 125%! 21.3333 6.3333 2.6667 1 16 ], 1e-4 ); 126 127%!test 128%! sz = 10; 129%! P = reshape( 1:sz^2, sz, sz ); 130%! normP = repmat(sum(P,2),1,columns(P)); 131%! P = P./normP; 132%! M = dtmcfpt(P); 133%! for i=1:rows(P) 134%! for j=1:columns(P) 135%! assert( M(i,j), 1 + dot(P(i,:), M(:,j)) - P(i,j)*M(j,j), 1e-8); 136%! endfor 137%! endfor 138 139## "Rat maze" problem (p. 453 of [GrSn97]); 140%!test 141%! P = zeros(9,9); 142%! P(1,[2 4]) = .5; 143%! P(2,[1 5 3]) = 1/3; 144%! P(3,[2 6]) = .5; 145%! P(4,[1 5 7]) = 1/3; 146%! P(5,[2 4 6 8]) = 1/4; 147%! P(6,[3 5 9]) = 1/3; 148%! P(7,[4 8]) = .5; 149%! P(8,[7 5 9]) = 1/3; 150%! P(9,[6 8]) = .5; 151%! M = dtmcfpt(P); 152%! assert( M(1:9 != 5,5)', [6 5 6 5 5 6 5 6], 100*eps ); 153 154