/* Macsyma commands to exponentiate the rate matrices A of examples 1
 * and 2 in Chapter 3, section 3.2, p.70 of the text.
 * AUTHOR: M.V.Wickerhauser
 * DATE: 2018-02-13
 * TEXT: Gregory F. Lawler, "Introduction to Stochastic Processes,"
 * second edition, ISBN 978-1-58488-651-8, Chapman and Hall/CRC, 2006.
 */

/* NOTE: the matrix product of A and B in Macsyma is A.B
 * while A*B gives the coordinate-wise product.
 */ 


load("eigen");			/* Make sure to have this extra Macsyma package */

A:matrix([-1,1],
	 [2,-2]);		/*** 2x2 rate matrix for Example 1 ***/

/*** DIAGONALIZTION ***/ 
/* Direct diagonalization of A using the eigenvectors of A: */ 
eigenvectors(A);		/* Compute eigenvalues and eigenvectors of A */
Q:transpose(matrix([1,1],[1,-2])); /* Matrix whose columns are eigenvectors of A */
D:invert(Q).A.Q;		/* Diagonalization of A gives diagonal matrix D */

/* Alternative diagonalization method available in Macsyma's "eigen" package: */ 
similaritytransform(A);		/* Compute eigenvalues and unit eigenvectors */
rightmatrix;			/* Matrix whose columns are unit eigenvectors */
leftmatrix;			/* Inverse of "rightmatrix" */
leftmatrix.A.rightmatrix;	/* Diagonalization of A. */

/*** EXPONENTIATION ***/ 
/* Exponentiate the matrix t*A for time t, and let P=exp(t*A): */ 
P:expand(matrixexp(t*A));	/* Separate the constant and exponential coefficients */

/* Find the stationary distribution from the long-time behavior of P 
 * by substituting 1000 for t in P, which should be long enough.  Then
 * evaluate the floating-point approximation so that the exponentially
 * decaying terms are truncated to zero. */
float(subst(1000,t,P));




A:matrix([-1,1,0,0],
	 [1,-3,1,1],
	 [0,1,-2,1],
	 [0,1,1,-2]);		/*** 4x4 rate matrix for Example 2 ***/

/*** DIAGONALIZTION ***/ 
/* Direct diagonalization of A using the eigenvectors of A: */ 
eigenvectors(A);		/* Compute eigenvalues and eigenvectors of A */

/* Alternative diagonalization method available in Macsyma's "eigen" package: */ 
similaritytransform(A);		/* Compute eigenvalues and unit eigenvectors */
rightmatrix;			/* Matrix whose columns are unit eigenvectors */
leftmatrix;			/* Inverse of "rightmatrix" */
leftmatrix.A.rightmatrix;	/* Diagonalization of A. */

/*** EXPONENTIATION ***/ 
/* Exponentiate the matrix t*A for time t, and let P=exp(t*A): */ 
P:expand(matrixexp(t*A));	/* Separate the constant and exponential coefficients */

/* Find the stationary distribution from the long-time behavior of P 
 * by substituting 1000 for t in P, which should be long enough.  Then
 * evaluate the floating-point approximation so that the exponentially
 * decaying terms are truncated to zero. */
float(subst(1000,t,P));


/*** WAITING TIMES ***/
/* ...to get to the last state from the other states */
Atilde: submatrix(4,A,4);	/* Remove last row 4 and last column 4 from A */
b: -invert(Atilde).matrix([1],[1],[1]); /* Equation from p.74 of the text. */