1## Copyright (C) 2012, 2014, 2016, 2018, 2019 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{t} @var{N} @var{B}] =} dtmcmtta (@var{P}) 21## @deftypefnx {Function File} {[@var{t} @var{N} @var{B}] =} dtmcmtta (@var{P}, @var{p0}) 22## 23## @cindex mean time to absorption, DTMC 24## @cindex absorption probabilities, DTMC 25## @cindex fundamental matrix 26## @cindex DTMC 27## @cindex discrete time Markov chain 28## @cindex Markov chain, discrete time 29## 30## Compute the expected number of steps before absorption for a 31## DTMC with state space @math{@{1, @dots{}, N@}} 32## and transition probability matrix @var{P}. 33## 34## @strong{INPUTS} 35## 36## @table @code 37## 38## @item @var{P}(i,j) 39## @math{N \times N} transition probability matrix. 40## @code{@var{P}(i,j)} is the transition probability from state 41## @math{i} to state @math{j}. 42## 43## @item @var{p0}(i) 44## Initial state occupancy probabilities (vector of length @math{N}). 45## 46## @end table 47## 48## @strong{OUTPUTS} 49## 50## @table @code 51## 52## @item @var{t} 53## @itemx @var{t}(i) 54## When called with a single argument, @var{t} is a vector of length 55## @math{N} such that @code{@var{t}(i)} is the expected number of steps 56## before being absorbed in any absorbing state, starting from state 57## @math{i}; if @math{i} is absorbing, @code{@var{t}(i) = 0}. When 58## called with two arguments, @var{t} is a scalar, and represents the 59## expected number of steps before absorption, starting from the initial 60## state occupancy probability @var{p0}. 61## 62## @item @var{N}(i) 63## @itemx @var{N}(i,j) 64## When called with a single argument, @var{N} is the @math{N \times N} 65## fundamental matrix for @var{P}. @code{@var{N}(i,j)} is the expected 66## number of visits to transient state @var{j} before absorption, if the 67## system started in transient state @var{i}. The initial state is counted 68## if @math{i = j}. When called with two arguments, @var{N} is a vector 69## of length @math{N} such that @code{@var{N}(j)} is the expected number 70## of visits to transient state @var{j} before absorption, given initial 71## state occupancy probability @var{P0}. 72## 73## @item @var{B}(i) 74## @itemx @var{B}(i,j) 75## When called with a single argument, @var{B} is a @math{N \times N} 76## matrix where @code{@var{B}(i,j)} is the probability of being 77## absorbed in state @math{j}, starting from transient state @math{i}; 78## if @math{j} is not absorbing, @code{@var{B}(i,j) = 0}; if @math{i} 79## is absorbing, @code{@var{B}(i,i) = 1} and @code{@var{B}(i,j) = 0} 80## for all @math{i \neq j}. When called with two arguments, @var{B} is 81## a vector of length @math{N} where @code{@var{B}(j)} is the 82## probability of being absorbed in state @var{j}, given initial state 83## occupancy probabilities @var{p0}. 84## 85## @end table 86## 87## @strong{REFERENCES} 88## 89## @itemize 90## @item Grinstead, Charles M.; Snell, J. Laurie (July 91## 1997). @cite{Introduction to Probability}, Ch. 11: Markov 92## Chains. American Mathematical Society. ISBN 978-0821807491. 93## @end itemize 94## 95## @seealso{ctmcmtta} 96## 97## @end deftypefn 98 99## Author: Moreno Marzolla <moreno.marzolla(at)unibo.it> 100## Web: http://www.moreno.marzolla.name/ 101 102function [t N B] = dtmcmtta( P, p0 ) 103 104 persistent epsilon = 10*eps; 105 106 if ( nargin < 1 || nargin > 2 ) 107 print_usage(); 108 endif 109 110 [K err] = dtmcchkP(P); 111 112 (K>0) || ... 113 error(err); 114 115 if ( nargin == 2 ) 116 ( isvector(p0) && length(p0) == K && all(p0>=0) && abs(sum(p0)-1.0)<epsilon ) || ... 117 error( "p0 must be a state occupancy probability vector" ); 118 endif 119 120 ## identify transient states 121 tr = find(diag(P) < 1); 122 ab = find(diag(P) == 1); 123 k = length(tr); # number of transient states 124 if ( k == K ) 125 error("There are no absorbing states"); 126 endif 127 128 N = B = zeros(size(P)); 129 t = zeros(1,rows(P)); 130 131 tmpN = inv(eye(k) - P(tr,tr)); # matrix N = (I-Q)^-1 132 N(tr,tr) = tmpN; 133 R = P(tr,ab); 134 135 res = tmpN * ones(k,1); 136 t(tr) = res; 137 138 tmp = tmpN*R; 139 B(tr,ab) = tmp; 140 ## set B(i,i) = 1 for all absorbing states i 141 dd = diag(B); 142 dd(ab) = 1; 143 B(1:K+1:end) = dd; 144 145 if ( nargin == 2 ) 146 t = dot(p0,t); 147 N = p0*N; 148 B = p0*B; 149 endif 150 151endfunction 152%!test 153%! fail( "dtmcmtta(1,2,3)" ); 154%! fail( "dtmcmtta()" ); 155 156%!test 157%! P = dtmcbd([0 .5 .5 .5], [.5 .5 .5 0]); 158%! [t N B] = dtmcmtta(P); 159%! assert( t, [0 3 4 3 0], 10*eps ); 160%! assert( B([2 3 4],[1 5]), [3/4 1/4; 1/2 1/2; 1/4 3/4], 10*eps ); 161%! assert( B(1,1), 1 ); 162%! assert( B(5,5), 1 ); 163 164%!test 165%! P = dtmcbd([0 .5 .5 .5], [.5 .5 .5 0]); 166%! [t N B] = dtmcmtta(P); 167%! assert( t(3), 4, 10*eps ); 168%! assert( B(3,1), 0.5, 10*eps ); 169%! assert( B(3,5), 0.5, 10*eps ); 170 171## Example on p. 422 of [GrSn97] 172%!test 173%! P = dtmcbd([0 .5 .5 .5 .5], [.5 .5 .5 .5 0]); 174%! [t N B] = dtmcmtta(P); 175%! assert( t(2:5), [4 6 6 4], 100*eps ); 176%! assert( B(2:5,1), [.8 .6 .4 .2]', 100*eps ); 177%! assert( B(2:5,6), [.2 .4 .6 .8]', 100*eps ); 178 179## Compute the probability of completing the "snakes and ladders" 180## game in n steps, for various values of n. Also, computes the expected 181## number of steps which are necessary to complete the game. 182## Source: 183## http://mapleta.emich.edu/aross15/coursepack3419/419/ch-04/chutes-and-ladders.pdf 184%!demo 185%! n = 6; 186%! P = zeros(101,101); 187%! for j=0:(100-n) 188%! i=1:n; 189%! P(1+j,1+j+i) = 1/n; 190%! endfor 191%! for j=(101-n):100 192%! P(1+j,1+j) = (n-100+j)/n; 193%! endfor 194%! for j=(101-n):100 195%! i=1:(100-j); 196%! P(1+j,1+j+i) = 1/n; 197%! endfor 198%! Pstar = P; 199%! ## setup snakes and ladders 200%! SL = [1 38; ... 201%! 4 14; ... 202%! 9 31; ... 203%! 16 6; ... 204%! 21 42; ... 205%! 28 84; ... 206%! 36 44; ... 207%! 47 26; ... 208%! 49 11; ... 209%! 51 67; ... 210%! 56 53; ... 211%! 62 19; ... 212%! 64 60; ... 213%! 71 91; ... 214%! 80 100; ... 215%! 87 24; ... 216%! 93 73; ... 217%! 95 75; ... 218%! 98 78 ]; 219%! for ii=1:rows(SL); 220%! i = SL(ii,1); 221%! j = SL(ii,2); 222%! Pstar(1+i,:) = 0; 223%! for k=0:100 224%! if ( k != i ) 225%! Pstar(1+k,1+j) = P(1+k,1+j) + P(1+k,1+i); 226%! endif 227%! endfor 228%! Pstar(:,1+i) = 0; 229%! endfor 230%! Pstar += diag( 1-sum(Pstar,2) ); 231%! # spy(Pstar); pause 232%! nsteps = 250; # number of steps 233%! Pfinish = zeros(1,nsteps); # Pfinish(i) = probability of finishing after step i 234%! pstart = zeros(1,101); pstart(1) = 1; pn = pstart; 235%! for i=1:nsteps 236%! pn = pn*Pstar; 237%! Pfinish(i) = pn(101); # state 101 is the ending (absorbing) state 238%! endfor 239%! f = dtmcmtta(Pstar,pstart); 240%! printf("Average number of steps to complete 'snakes and ladders': %f\n", f ); 241%! plot(Pfinish,"linewidth",2); 242%! line([f,f],[0,1]); 243%! text(f*1.1,0.2,["Mean Time to Absorption (" num2str(f) ")"]); 244%! xlabel("Step number (n)"); 245%! title("Probability of finishing 'snakes and ladders' before step n"); 246 247## "Rat maze" problem (p. 453 of [GrSn97]); 248%!test 249%! P = zeros(9,9); 250%! P(1,[2 4]) = .5; 251%! P(2,[1 5 3]) = 1/3; 252%! P(3,[2 6]) = .5; 253%! P(4,[1 5 7]) = 1/3; 254%! P(5,:) = 0; P(5,5) = 1; 255%! P(6,[3 5 9]) = 1/3; 256%! P(7,[4 8]) = .5; 257%! P(8,[7 5 9]) = 1/3; 258%! P(9,[6 8]) = .5; 259%! t = dtmcmtta(P); 260%! assert( t, [6 5 6 5 0 5 6 5 6], 10*eps ); 261 262