Re: General 3-state stochastic matrix
- To: mathgroup at smc.vnet.net
- Subject: [mg58227] Re: General 3-state stochastic matrix
- From: Maxim <ab_def at prontomail.com>
- Date: Thu, 23 Jun 2005 05:34:17 -0400 (EDT)
- References: <d98p16$eeo$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
On Tue, 21 Jun 2005 10:09:42 +0000 (UTC), Virgil Stokes
<virgil.stokes at it.uu.se> wrote:
> I have tried to find the limit (as n, the power of the matrix, goes to
> infinity) for the general 3-state stochastic matrix using the following
> code:
>
...
>
> However, it does not find a symbolic solution. I would appreciate it
> greatly if someone else could look at this and see if they are able to
> get a symbolic solution. Warning! this can take considerable CPU time.
>
> Note, for a general 2-state stochastic matrix, the above approach works
> fine.
>
> --Thanks, V. Stokes
>
Actually, T is a transpose of a stochastic matrix. Mathematica will most
likely have problems with this limit, because even in cases like
Limit[(a/2)^n, n -> Infinity, Assumptions -> a > 2]
Limit is just left unevaluated. One method is to use the properties of X =
T^Infinity:
In[1]:=
T = {{1 - a2 - a3, a2, a3}, {b1, 1 - b1 - b3, b3}, {c1, c2, 1 - c1 - c2}};
X = Array[x, {3, 3}];
In[3]:=
sol = First@ Solve[T.X == X.T == X, Variables@ X];
sol2 = First@ Solve[Eigenvalues@ X == {0, 0, 1} /. sol, Variables@ X];
X /. sol /. sol2 // Simplify
Out[5]=
{{(b3*c1 + b1*(c1 + c2))/(b1*c1 + b3*c1 + b1*c2 + a3*(b1 + b3 + c2)
+ a2*(b3 + c1 + c2)), (a3*c2 + a2*(c1 + c2))/(b1*c1 + b3*c1 + b1*c2
+ a3*(b1 + b3 + c2) + a2*(b3 + c1 + c2)), (a2*b3 + a3*(b1 + b3))/(b1*c1
+ b3*c1 + b1*c2 + a3*(b1 + b3 + c2) + a2*(b3 + c1 + c2))}, {(b3*c1
+ b1*(c1 + c2))/(b1*c1 + b3*c1 + b1*c2 + a3*(b1 + b3 + c2) + a2*(b3 + c1
+ c2)), (a3*c2 + a2*(c1 + c2))/(b1*c1 + b3*c1 + b1*c2 + a3*(b1 + b3 + c2)
+ a2*(b3 + c1 + c2)), (a2*b3 + a3*(b1 + b3))/(b1*c1 + b3*c1 + b1*c2
+ a3*(b1 + b3 + c2) + a2*(b3 + c1 + c2))}, {(b3*c1 + b1*(c1 + c2))/(b1*c1
+ b3*c1 + b1*c2 + a3*(b1 + b3 + c2) + a2*(b3 + c1 + c2)), (a3*c2 + a2*(c1
+ c2))/(b1*c1 + b3*c1 + b1*c2 + a3*(b1 + b3 + c2) + a2*(b3 + c1 + c2)),
(a2*b3 + a3*(b1 + b3))/(b1*c1 + b3*c1 + b1*c2 + a3*(b1 + b3 + c2) + a2*(b3
+ c1 + c2))}}
Or
In[6]:=
Module[{s, j},
{s, j} = JordanDecomposition[T];
s.Replace[j, {1 -> 1, _ -> 0}, {2}].Inverse[s] == %5 // Simplify
]
Out[6]=
True
I think this will hold for degenerate cases too, as an appropriate limit.
Maxim Rytin
m.r at inbox.ru