changeset 9733:1a4d4ff140c7 octave-forge

bug fixes
author mmarzolla
date Fri, 16 Mar 2012 14:08:44 +0000
parents 0f2ce9bc2525
children d181b0e4c274
files main/queueing/inst/dtmc.m main/queueing/inst/dtmc_check_P.m main/queueing/inst/dtmc_fpt.m main/queueing/inst/dtmc_is_irreducible.m main/queueing/inst/dtmc_mtta.m
diffstat 5 files changed, 60 insertions(+), 32 deletions(-) [+]
line wrap: on
line diff
--- a/main/queueing/inst/dtmc.m	Fri Mar 16 11:24:04 2012 +0000
+++ b/main/queueing/inst/dtmc.m	Fri Mar 16 14:08:44 2012 +0000
@@ -24,14 +24,16 @@
 ## @cindex Discrete time Markov chain
 ## @cindex Markov chain, stationary probabilities
 ## @cindex Stationary probabilities
+## @cindex Markov chain, transient probabilities
+## @cindex Transient probabilities
 ##
-## With a single argument, compute the steady-state probability vector
-## @code{@var{p}(1), @dots{}, @var{p}(N)} for a
-## Discrete-Time Markov Chain given the @math{N \times N} transition
-## probability matrix @var{P}. With three arguments, compute the
-## probability vector @code{@var{p}(1), @dots{}, @var{p}(N)}
-## after @var{n} steps, given initial probability vector @var{p0} at
-## time 0.
+## Compute steady-state or transient state occupancy probabilities for a
+## Discrete-Time Markov Chain. With a single argument, compute the
+## steady-state occupancy probability vector @code{@var{p}(1), @dots{},
+## @var{p}(N)} given the @math{N \times N} transition probability matrix
+## @var{P}. With three arguments, compute the state occupancy
+## probabilities @code{@var{p}(1), @dots{}, @var{p}(N)} after @var{n}
+## steps, given initial occupancy probability vector @var{p0}.
 ##
 ## @strong{INPUTS}
 ##
@@ -145,6 +147,15 @@
 %! p = dtmc(P, 100, [1 0]);
 %! assert( plim, p, 1e-5 );
 
+%!test
+%! P = [0 1 0 0 0; \
+%!      .25 0 .75 0 0; \
+%!      0 .5 0 .5 0; \
+%!      0 0 .75 0 .25; \
+%!      0 0 0 1 0 ];
+%! p = dtmc(P);
+%! assert( p, [.0625 .25 .375 .25 .0625], 10*eps );
+
 %!demo
 %! a = 0.2;
 %! b = 0.15;
--- a/main/queueing/inst/dtmc_check_P.m	Fri Mar 16 11:24:04 2012 +0000
+++ b/main/queueing/inst/dtmc_check_P.m	Fri Mar 16 14:08:44 2012 +0000
@@ -17,14 +17,14 @@
 
 ## -*- texinfo -*-
 ##
-## @deftypefn {Function File} {[@var{result} @var{err}] =} dtmc_check_P (@var{P})
+## @deftypefn {Function File} {[@var{r} @var{err}] =} dtmc_check_P (@var{P})
 ##
 ## @cindex Markov chain, discrete time
 ##
-## If @var{P} is a valid transition probability matrix, return
-## the size (number of rows or columns) of @var{P}. If @var{P} is not
-## a transition probability matrix, set @var{result} to zero, and
-## @var{err} to an appropriate error string.
+## Check if @var{P} is a valid transition probability matrix. If @var{P}
+## is valid, @var{r} is the size (number of rows or columns) of @var{P}.
+## If @var{P} is not a transition probability matrix, @var{r} is set to
+## zero, and @var{err} to an appropriate error string.
 ##
 ## @end deftypefn
 
@@ -40,19 +40,19 @@
   endif
 
   result = 0;
+  err = "";
 
   if ( !issquare(P) )
     err = "P is not a square matrix";
     return;
   endif
   
-  if (  any(any(P <0)) || norm( sum(P,2) - 1, "inf" ) > epsilon )
+  if (  any(any(P<-epsilon)) || norm( sum(P,2) - 1, "inf" ) > epsilon )
     err = "P is not a stochastic matrix";
     return;
   endif
 
   result = rows(P);
-  err = "";
 endfunction
 %!test
 %! [r err] = dtmc_check_P( [1 1 1; 1 1 1] );
--- a/main/queueing/inst/dtmc_fpt.m	Fri Mar 16 11:24:04 2012 +0000
+++ b/main/queueing/inst/dtmc_fpt.m	Fri Mar 16 14:08:44 2012 +0000
@@ -84,6 +84,10 @@
   ( N>0 ) || \
       error(err);
 
+  if ( any(diag(P) == 1) ) 
+    error("Cannot compute first passage times for absorbing chains");
+  endif
+
   if ( nargin == 1 )   
     M = zeros(N,N);
     ## M(i,j) = 1 + sum_{k \neq j} P(i,k) M(k,j)
@@ -119,6 +123,13 @@
 %!       0.1 0.0 0.9; \
 %!       0.9 0.1 0.0 ];
 %! M = dtmc_fpt(P);
+%! w = dtmc(P);
+%! N = rows(P);
+%! W = repmat(w,N,1);
+%! Z = inv(eye(N)-P+W);
+%! M1 = (repmat(diag(Z)',1,N) - Z) ./ repmat(w',1,N);
+%! assert(M, M1);
+
 
 %!shared P
 %! P = [ 0.0 0.9 0.1; \
@@ -138,10 +149,9 @@
 %!test
 %! m = dtmc_fpt(P, 1, [2 3]);
 
-## FIXME: check this (matrix not ergodic???)
-%!xtest
+%!test
 %! P = dtmc_bd([1 1 1], [ 0 0 0] );
-%! dtmc_fpt(P);
+%! fail( "dtmc_fpt(P)", "absorbing" );
 
 ## Example on p. 461 of
 ## http://www.cs.virginia.edu/~gfx/Courses/2006/DataDriven/bib/texsyn/Chapter11.pdf
--- a/main/queueing/inst/dtmc_is_irreducible.m	Fri Mar 16 11:24:04 2012 +0000
+++ b/main/queueing/inst/dtmc_is_irreducible.m	Fri Mar 16 14:08:44 2012 +0000
@@ -21,10 +21,11 @@
 ##
 ## @cindex Markov chain, discrete time
 ## @cindex Discrete time Markov chain
+## @cindex Irreducible Markov chain
 ##
-## Check if @var{P} is irreducible, and identify strongly connected
-## components in the transition graph of the discrete-time Markov chain
-## with transition probability matrix @var{P}.
+## Check if @var{P} is irreducible, and identify Strongly Connected
+## Components (SCC) in the transition graph of the DTMC with transition
+## probability matrix @var{P}.
 ##
 ## @strong{INPUTS}
 ##
@@ -45,10 +46,10 @@
 ## 1 if @var{P} irreducible, 0 otherwise.
 ##
 ## @item s
-## @code{@var{s}(i)} is the strongly connected component that state @math{i}
-## belongs to. Strongly connected components are numbered as 1, 2,
-## @dots{}. If the graph is strongly connected, then there is a single
-## SCC and the predicate @code{all(s == 1)} evaluates to true.
+## @code{@var{s}(i)} is the SCC that state @math{i} belongs to. SCCs are
+## numbered as 1, 2, @dots{}. If the graph is strongly connected, then
+## there is a single SCC and the predicate @code{all(s == 1)} evaluates
+## to true.
 ##
 ## @end table
 ##
@@ -59,19 +60,23 @@
 
 function [r s] = dtmc_is_irreducible( P )
 
-  persistent epsilon = 10*eps;
-
   if ( nargin != 1 )
     print_usage();
   endif
 
-  # dtmc_check_P(P);
-
+  [N err] = dtmc_check_P(P);
+  if ( N == 0 ) 
+    error(err);
+  endif
   s = __scc(P);
   r = (max(s) == 1);
 
 endfunction
 %!test
+%! P = [0 .5 0; 0 0 0];
+%! fail( "dtmc_is_irresudible(P)" );
+
+%!test
 %! P = [0 1 0; 0 .5 .5; 0 1 0];
 %! [r s] = dtmc_is_irreducible(P);
 %! assert( r == 0 );
--- a/main/queueing/inst/dtmc_mtta.m	Fri Mar 16 11:24:04 2012 +0000
+++ b/main/queueing/inst/dtmc_mtta.m	Fri Mar 16 14:08:44 2012 +0000
@@ -22,6 +22,7 @@
 ##
 ## @cindex Markov chain, disctete time
 ## @cindex Mean time to absorption
+## @cindex Absorption probabilities
 ##
 ## Compute the expected number of steps before absorption for the 
 ## DTMC with @math{N \times N} transition probability matrix @var{P}.
@@ -53,9 +54,10 @@
 ## @item B
 ## When called with a single argument, @var{B} is a @math{N \times N}
 ## matrix where @code{@var{B}(i,j)} is the probability of being absorbed
-## in state @math{j}, starting from state @math{i}; if @math{j} is not
-## absorbing, @code{@var{B}(i,j) = 0}; if @math{i} is absorbing, then
-## @code{@var{B}(i,i) = 1}. When called with two arguments, @var{B} is a
+## in state @math{j}, starting from transient state @math{i}; if
+## @math{j} is not absorbing, @code{@var{B}(i,j) = 0}; if @math{i} is
+## absorbing, then @code{@var{B}(i,i) = 1}.
+## When called with two arguments, @var{B} is a
 ## vector with @math{N} elements where @code{@var{B}(j)} is the
 ## probability of being absorbed in state @var{j}, given initial state
 ## occupancy probabilities @var{p0}.
@@ -82,7 +84,7 @@
   
   if ( nargin == 2 )
     ( isvector(p0) && length(p0) == K && all(p0>=0) && abs(sum(p0)-1.0)<epsilon ) || \
-	usage( "p0 must be a probability vector" );
+	usage( "p0 must be a state occupancy probability vector" );
   endif