-— Function File: L = ctmc_exps (Q, tt, p )
+— Function File: L = ctmc_exps (Q, t, p )
— Function File: L = ctmc_exps (Q, p)

-With three arguments, compute the expected time L(t,j) -spent in each state j during the time interval -[0,tt(t)), assuming that at time 0 the state occupancy -probability was p. With two arguments, compute the expected -time L(j) spent in each state j until absorption. +With three arguments, compute the expected times L(i) +spent in each state i during the time interval +[0,t], assuming that the state occupancy probabilities +at time 0 are p. With two arguments, compute the expected time +L(i) spent in each state i until absorption.
INPUTS @@ -1118,10 +1118,7 @@ ≤ i \neq j ≤ N. The matrix Q must also satisfy the condition \sum_j=1^N Q_ij = 0. -
tt
This parameter is a vector used for numerical integration. The first -element tt(1) must be 0, and the last element -tt(end) must be the upper bound of the interval -[0,t) of interest (tt(end) == t). +
t
Time
p
Initial occupancy probability vector; p(i) is the probability the system is in state i at time 0, i = 1, @@ -1132,14 +1129,12 @@
OUTPUTS
-
L
If this function is called with three arguments, L is a matrix -of size [length(tt), N] where L(t,j) is the -expected time spent in state j during the interval -[0,tt(t)]. If this function is called with two -arguments, L is a vector with N elements where -L(j) is the expected time spent in state j until -absorption, if j is a transient state. If j -is an absorbing state, L(j) = 0. +
L
If this function is called with three arguments, L(i) is +the expected time spent in state j during the interval +[0,t]. If this function is called with two arguments +L(i) is the expected time spent in state i until +absorption (if i is a transient state), or zero +(if i is an absorbing state).
@@ -1185,13 +1180,13 @@

-— Function File: M = ctmc_taexps (Q, tt, p)
+— Function File: M = ctmc_taexps (Q, t, p)

-Compute the time-averaged sojourn time M(t,j), -defined as the fraction of the time interval [0,tt(t)) spent in -state j, assuming that at time 0 the state occupancy -probability was p. +Compute the time-averaged sojourn time M(i), +defined as the fraction of the time interval [0,t] spent in +state i, assuming that the state occupancy probabilities at +time 0 are p.
INPUTS @@ -1199,15 +1194,9 @@
Q
Infinitesimal generator matrix. Q(i,j) is the transition rate from state i to state j, 1 ≤ i \neq j ≤ N. The -matrix Q must also satisfy the condition sum(Q,2) == 0 - -
tt
This parameter is a vector used for numerical integration of the -sujourn time. The first element tt(1) must be slightly -larger than 0, and the -last element tt(end) must be the upper limit of the -interval [0,t) of interest (tt(end) == t). -This vector is used by the ODE solver to compute the solution -M. +matrix Q must also satisfy the condition \sum_j=1^N Q_ij = 0 + +
t
Time
p
p(i) is the probability that, at time 0, the system was in state i, for all i = 1, ..., N @@ -1217,10 +1206,9 @@
OUTPUTS
-
M
M(t,j) is the expected fraction of time spent in state -j during the interval [0,tt(t)) assuming that the state -occupancy probability at time zero was p. 1 ≤ -t ≤ length(tt) +
M
M(i) is the expected fraction of time spent in state +i during the interval [0,t] assuming that the state +occupancy probability at time zero is p.
@@ -1279,7 +1267,7 @@
Compute the Mean-Time to Absorption (MTTA) of the CTMC described by the infinitesimal generator matrix Q, starting from initial -occupancy probability p. If there are no absorbing states, this +occupancy probabilities p. If there are no absorbing states, this function fails with an error.
INPUTS diff -r dbea4e0ca957 -r c2e1990614a2 main/queueing/doc/queueing.pdf Binary file main/queueing/doc/queueing.pdf has changed diff -r dbea4e0ca957 -r c2e1990614a2 main/queueing/inst/ctmc_exps.m --- a/main/queueing/inst/ctmc_exps.m Sat Mar 10 19:46:27 2012 +0000 +++ b/main/queueing/inst/ctmc_exps.m Sat Mar 10 20:27:25 2012 +0000 @@ -17,17 +17,17 @@ ## -*- texinfo -*- ## -## @deftypefn {Function File} {@var{L} =} ctmc_exps (@var{Q}, @var{tt}, @var{p} ) +## @deftypefn {Function File} {@var{L} =} ctmc_exps (@var{Q}, @var{t}, @var{p} ) ## @deftypefnx {Function File} {@var{L} =} ctmc_exps (@var{Q}, @var{p}) ## ## @cindex Markov chain, continuous time ## @cindex Expected sojourn time ## -## 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. +## With three arguments, compute the expected times @code{@var{L}(i)} +## spent in each state @math{i} during the time interval +## @math{[0,t]}, assuming that the state occupancy probabilities +## at time 0 are @var{p}. With two arguments, compute the expected time +## @code{@var{L}(i)} spent in each state @math{i} until absorption. ## ## @strong{INPUTS} ## @@ -39,11 +39,8 @@ ## @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 -## element @code{@var{tt}(1)} must be 0, and the last element -## @code{@var{tt}(end)} must be the upper bound of the interval -## @math{[0,t)} of interest (@code{@var{tt}(end) == @math{t}}). +## @item t +## Time ## ## @item p ## Initial occupancy probability vector; @code{@var{p}(i)} is the @@ -57,14 +54,12 @@ ## @table @var ## ## @item L -## 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}. +## If this function is called with three arguments, @code{@var{L}(i)} is +## the expected time spent in state @math{j} during the interval +## @math{[0,t]}. If this function is called with two arguments +## @code{@var{L}(i)} is the expected time spent in state @math{i} until +## absorption (if @math{i} is a transient state), or zero +## (if @var{i} is an absorbing state). ## ## @end table ## @@ -85,12 +80,12 @@ usage( "Q must be a square matrix" ); ( norm( sum(Q,2), "inf" ) < epsilon ) || \ - error( "Q is not an infinitesimal generator matrix" ); + usage( "Q is not an infinitesimal generator matrix" ); if ( nargin == 2 ) p = varargin{1}; else - tt = varargin{1}; + t = varargin{1}; p = varargin{2}; endif @@ -98,15 +93,27 @@ usage( "p must be a probability vector" ); 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 ); + if ( isscalar(t) ) + (t >= 0 ) || \ + usage( "t must be >= 0" ); + ## F(x) are the transient state occupancy probabilities at time x. + ## It is known that F(x) = p*expm(Q*x) (see function ctmc()). + F = @(x) (p*expm(Q*x)); + L = quadv(F,0,t); + else + ## FIXME: deprecate this? + ( isvector(t) && abs(t(1)) < epsilon ) || \ + usage( "t must be a vector, and t(1) must be 0.0" ); + t = t(:)'; # 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)), t ); + endif else #{ + ## This code is left for information only + ## 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 @@ -138,6 +145,10 @@ L(nzrows) = LN; endif endfunction +%!test +%! Q = [-1 1; 1 -1]; +%! L = ctmc_exps(Q,10,[1 0]); +%! L = ctmc_exps(Q,linspace(0,10,100),[1 0]); %!demo %! lambda = 0.5; @@ -146,18 +157,16 @@ %! death = zeros(1,N-1); %! Q = diag(birth,1)+diag(death,-1); %! Q -= diag(sum(Q,2)); -%! tt = linspace(0,10,100); +%! t = 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); +%! L = zeros(length(t),N); +%! for i=1:length(t) +%! L(i,:) = ctmc_exps(Q,t(i),p0); +%! endfor +%! plot( t, L(:,1), ";State 1;", "linewidth", 2, \ +%! t, L(:,2), ";State 2;", "linewidth", 2, \ +%! t, L(:,3), ";State 3;", "linewidth", 2, \ +%! t, L(:,4), ";State 4;", "linewidth", 2 ); %! legend("location","northwest"); %! xlabel("Time"); %! ylabel("Expected sojourn time"); diff -r dbea4e0ca957 -r c2e1990614a2 main/queueing/inst/ctmc_mtta.m --- a/main/queueing/inst/ctmc_mtta.m Sat Mar 10 19:46:27 2012 +0000 +++ b/main/queueing/inst/ctmc_mtta.m Sat Mar 10 20:27:25 2012 +0000 @@ -24,7 +24,7 @@ ## ## 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 +## occupancy probabilities @var{p}. If there are no absorbing states, this ## function fails with an error. ## ## @strong{INPUTS} @@ -71,8 +71,8 @@ N = rows(Q); - all( abs( sum(Q,2) ) < epsilon ) || \ - usage( "Q is not an infinitesimal generator matrix" ); + ( norm( sum(Q,2), "inf" ) < epsilon ) || \ + usage( "Q must be an infinitesimal generator matrix" ); ( isvector(p) && length(p) == N && all(p>=0) && abs(sum(p)-1.0)=0) && abs(sum(p)-1.0)= 0) || \ + usage( "t must be >= 0" ); + F = @(x) (p*expm(Q*x)); + M = quadv(F,0,t) / t; + else + ## FIXME: deprecate this? + t = t(:)'; # make t a row vector + p = p(:)'; # make p a row vector + ff = @(x,t) (((x')*(Q-eye(N)/t).+p/t)'); + fj = @(x,t) (Q-eye(N)/t); + M = lsode( {ff, fj}, zeros(size(p)), t ); + endif endfunction %!demo @@ -101,9 +102,12 @@ %! death = zeros(1,N-1); %! Q = diag(birth,1)+diag(death,-1); %! Q -= diag(sum(Q,2)); -%! t = linspace(1e-3,50,500); +%! t = linspace(1e-5,30,100); %! p = zeros(1,N); p(1)=1; -%! M = ctmc_taexps(Q,t,p); +%! M = zeros(length(t),N); +%! for i=1:length(t) +%! M(i,:) = ctmc_taexps(Q,t(i),p); +%! endfor %! plot(t, M(:,1), ";State 1;", "linewidth", 2, \ %! t, M(:,2), ";State 2;", "linewidth", 2, \ %! t, M(:,3), ";State 3;", "linewidth", 2, \