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