From: <mma...@us...> - 2012-03-10 16:00:52
|
Revision: 9807 http://octave.svn.sourceforge.net/octave/?rev=9807&view=rev Author: mmarzolla Date: 2012-03-10 16:00:45 +0000 (Sat, 10 Mar 2012) Log Message: ----------- Bug fix and enhancements in ctmc_exps() Modified Paths: -------------- trunk/octave-forge/main/queueing/ChangeLog trunk/octave-forge/main/queueing/NEWS trunk/octave-forge/main/queueing/inst/ctmc_exps.m trunk/octave-forge/main/queueing/inst/ctmc_mtta.m Modified: trunk/octave-forge/main/queueing/ChangeLog =================================================================== --- trunk/octave-forge/main/queueing/ChangeLog 2012-03-10 12:44:08 UTC (rev 9806) +++ trunk/octave-forge/main/queueing/ChangeLog 2012-03-10 16:00:45 UTC (rev 9807) @@ -1,3 +1,13 @@ +2012-02-XX Moreno Marzolla <mar...@cs...> + + * Version 1.0.X released + * Fixed bug in qnvisits() which made the function behave incorrectly + under particular degenerate cases. + * Fixed bug in ctmc_exps() (wrong initial value in call to lsode) + * Function ctmc_exps() can now also compute the expected sojourn time + until absorption for absorbing CTMCs. + * Miscellaneous fixes/improvements to the documentation + 2012-02-04 Moreno Marzolla <mar...@cs...> * Version 1.0.0 released under the name "queueing" (initial Modified: trunk/octave-forge/main/queueing/NEWS =================================================================== --- trunk/octave-forge/main/queueing/NEWS 2012-03-10 12:44:08 UTC (rev 9806) +++ trunk/octave-forge/main/queueing/NEWS 2012-03-10 16:00:45 UTC (rev 9807) @@ -1,3 +1,8 @@ +Summary of important user-visible changes for queueing-1.0.X + +** Function ctmc_exps() can now compute the expected sojourn time + until absorption for absorming CTMC + Summary of important user-visible changes for queueing-1.0.0 ------------------------------------------------------------------------------ Modified: trunk/octave-forge/main/queueing/inst/ctmc_exps.m =================================================================== --- trunk/octave-forge/main/queueing/inst/ctmc_exps.m 2012-03-10 12:44:08 UTC (rev 9806) +++ trunk/octave-forge/main/queueing/inst/ctmc_exps.m 2012-03-10 16:00:45 UTC (rev 9807) @@ -17,24 +17,27 @@ ## -*- texinfo -*- ## -## @deftypefn {Function File} {@var{L} =} ctmc_exps (@var{Q}, @var{tt}, @var{p}) +## @deftypefn {Function File} {@var{L} =} ctmc_exps (@var{Q}, @var{tt}, @var{p} ) +## @deftypefnx {Function File} {@var{L} =} ctmc_exps (@var{Q}, @var{p}) ## ## @cindex Markov chain, continuous time ## @cindex Expected sojourn time ## -## Compute the expected total time @code{@var{L}(t,j)} spent in state -## @math{j} during the time interval @code{[0,@var{tt}(t))}, assuming -## that at time 0 the state occupancy probability was @var{p}. +## With three arguments, compute the expected time @code{@var{L}(t,j)} +## spent in each state @math{j} during the time interval +## @code{[0,@var{tt}(t))}, assuming that at time 0 the state occupancy +## probability was @var{p}. With two arguments, compute the expected +## time @code{@var{L}(j}} spent in each state @math{j} until absorption. ## ## @strong{INPUTS} ## ## @table @var ## ## @item Q -## Infinitesimal generator matrix. @code{@var{Q}(i,j)} is the transition -## rate from state @math{i} to state @math{j}, -## @math{1 @leq{} i \neq j @leq{} N}. The matrix @var{Q} must also satisfy the -## condition @code{sum(@var{Q},2) == 0} +## @math{N \times N} infinitesimal generator matrix. @code{@var{Q}(i,j)} +## is the transition rate from state @math{i} to state @math{j}, @math{1 +## @leq{} i \neq j @leq{} N}. The matrix @var{Q} must also satisfy the +## condition @math{\sum_{j=1}^N Q_{ij} = 0}. ## ## @item tt ## This parameter is a vector used for numerical integration. The first @@ -43,8 +46,9 @@ ## @math{[0,t)} of interest (@code{@var{tt}(end) == @math{t}}). ## ## @item p -## @code{@var{p}(i)} is the probability that at time 0 the system was in -## state @math{i}, for all @math{i = 1, @dots{}, N} +## Initial occupancy probability vector; @code{@var{p}(i)} is the +## probability the system is in state @math{i} at time 0, @math{i = 1, +## @dots{}, N} ## ## @end table ## @@ -53,8 +57,14 @@ ## @table @var ## ## @item L -## @code{@var{L}(t,j)} is the expected time spent in state @math{j} -## during the interval @code{[0,@var{tt}(t))}. @code{1 @leq{} @var{t} @leq{} length(@var{tt})} +## If this function is called with three arguments, @var{L} is a matrix +## of size @code{[length(@var{tt}), N]} where @code{@var{L}(t,j)} is the +## expected time spent in state @math{j} during the interval +## @code{[0,@var{tt}(t)]}. If this function is called with two +## arguments, @var{L} is a vector with @math{N} elements where +## @code{@var{L}(j)} is the expected time spent in state @math{j} until +## absorption, if @math{j} is a transient state. If @math{j} +## is an absorbing state, @code{@var{L}(j) = 0}. ## ## @end table ## @@ -63,11 +73,11 @@ ## Author: Moreno Marzolla <marzolla(at)cs.unibo.it> ## Web: http://www.moreno.marzolla.name/ -function L = ctmc_exps( Q, tt, p ) +function L = ctmc_exps( Q, varargin ) persistent epsilon = 10*eps; - if ( nargin != 3 ) + if ( nargin < 2 || nargin > 3 ) print_usage(); endif @@ -77,16 +87,56 @@ ( norm( sum(Q,2), "inf" ) < epsilon ) || \ error( "Q is not an infinitesimal generator matrix" ); + if ( nargin == 2 ) + p = varargin{1}; + else + tt = varargin{1}; + p = varargin{2}; + endif + ( isvector(p) && length(p) == size(Q,1) && all(p>=0) && abs(sum(p)-1.0)<epsilon ) || \ usage( "p must be a probability vector" ); - ( isvector(tt) && abs(tt(1)) < epsilon ) || \ - usage( "tt must be a vector, and tt(1) must be 0.0" ); - tt = tt(:)'; # make tt a row vector - p = p(:)'; # make p a row vector - ff = @(x,t) (x(:)'*Q+p); - fj = @(x,t) (Q); - L = lsode( {ff, fj}, p, tt ); + if ( nargin == 3 ) + ( isvector(tt) && abs(tt(1)) < epsilon ) || \ + usage( "tt must be a vector, and tt(1) must be 0.0" ); + tt = tt(:)'; # make tt a row vector + p = p(:)'; # make p a row vector + ff = @(x,t) (x(:)'*Q+p); + fj = @(x,t) (Q); + L = lsode( {ff, fj}, zeros(size(p)), tt ); + else +#{ + ## F(t) are the transient state occupancy probabilities at time t. + ## It is known that F(t) = p*expm(Q*t) (see function ctmc()). + ## The expected times spent in each state until absorption can + ## be computed as the integral of F(t) from t=0 to t=inf + F = @(t) (p*expm(Q*t)); ## FIXME: this must be restricted to transient states ONLY!!!! + + ## Since function quadv does not support infinite integration + ## limits, we define a new function G(u) = F(tan(pi/2*u)) such that + ## the integral of G(u) on [0,1] is the integral of F(t) on [0, + ## +inf]. + G = @(u) (F(tan(pi/2*u))*pi/2*(1+tan(pi/2*u)**2)); + + L = quadv(G,0,1); +#} + ## Find nonzero rows. Nonzero rows correspond to transient states, + ## while zero rows are absorbing states. If there are no zero rows, + ## then the Markov chain does not contain absorbing states and we + ## raise an error + N = rows(Q); + nzrows = find( any( abs(Q) > epsilon, 2 ) ); + if ( length( nzrows ) == N ) + error( "There are no absorbing states" ); + endif + + QN = Q(nzrows,nzrows); + pN = p(nzrows); + LN = -pN*inv(QN); + L = zeros(1,N); + L(nzrows) = LN; + endif endfunction %!demo @@ -99,10 +149,25 @@ %! tt = linspace(0,10,100); %! p0 = zeros(1,N); p0(1)=1; %! L = ctmc_exps(Q,tt,p0); +%! #L2 = 0*L; +%! #for i=1:length(tt) +%! # L2(i,:) = quadv( @(t) (p0*expm(Q*t)) , 0, tt(i) ); +%! #endfor %! plot( tt, L(:,1), ";State 1;", "linewidth", 2, \ +%! # tt, L2(:,1), "+;State 1 (quadv);", "markersize", 8, \ %! tt, L(:,2), ";State 2;", "linewidth", 2, \ %! tt, L(:,3), ";State 3;", "linewidth", 2, \ %! tt, L(:,4), ";State 4 (absorbing);", "linewidth", 2); %! legend("location","northwest"); %! xlabel("Time"); %! ylabel("Expected sojourn time"); + +%!demo +%! lambda = 0.5; +%! N = 4; +%! birth = lambda*linspace(1,N-1,N-1); +%! death = 0*birth; +%! Q = diag(birth,1)+diag(death,-1); +%! Q -= diag(sum(Q,2)); +%! p0 = zeros(1,N); p0(1)=1; +%! L = ctmc_exps(Q,p0) Modified: trunk/octave-forge/main/queueing/inst/ctmc_mtta.m =================================================================== --- trunk/octave-forge/main/queueing/inst/ctmc_mtta.m 2012-03-10 12:44:08 UTC (rev 9806) +++ trunk/octave-forge/main/queueing/inst/ctmc_mtta.m 2012-03-10 16:00:45 UTC (rev 9807) @@ -22,9 +22,10 @@ ## @cindex Markov chain, continuous time ## @cindex Mean time to absorption ## -## Compute the Mean-Time to Absorption (MTTA) starting from initial -## occupancy probability @var{p} at time 0. If there are no absorbing -## states, this function fails with an error. +## Compute the Mean-Time to Absorption (MTTA) of the CTMC described by +## the infinitesimal generator matrix @var{Q}, starting from initial +## occupancy probability @var{p}. If there are no absorbing states, this +## function fails with an error. ## ## @strong{INPUTS} ## @@ -76,18 +77,7 @@ ( isvector(p) && length(p) == N && all(p>=0) && abs(sum(p)-1.0)<epsilon ) || \ usage( "p must be a probability vector" ); - ## Find nonzero rows. Nonzero rows correspond to transient states, - ## while zero rows are absorbing states. If there are no zero rows, - ## then the Markov chain does not contain absorbing states and we - ## raise an error - nzrows = find( any( abs(Q) > epsilon, 2 ) ); - if ( length( nzrows ) == N ) - error( "There are no absorbing states" ); - endif - - QN = Q(nzrows,nzrows); - pN = p(nzrows); - L = -pN*inv(QN); + L = ctmc_exps(Q,p); t = sum(L); endfunction %!test @@ -122,7 +112,7 @@ %! death = [ 3 4 5 ] * mu; %! Q = diag(death,-1); %! Q -= diag(sum(Q,2)); -%! t = ctmc_mtta(Q,[0 0 0 1]) +%! [t L] = ctmc_mtta(Q,[0 0 0 1]) %!demo %! N = 100; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |