Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2004
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2004

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: any tool/package to solve markov chain?

  • To: mathgroup at smc.vnet.net
  • Subject: [mg51948] Re: [mg51869] any tool/package to solve markov chain?
  • From: Daniel Lichtblau <danl at wolfram.com>
  • Date: Fri, 5 Nov 2004 02:18:59 -0500 (EST)
  • References: <200411040649.BAA18023@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

Yan ZHANG wrote:
> Now, I got a Markov chain and expressed as x=xP where x is the steady state
> probability vector and P is the transition matrix. Since the dimension of P
> is very big (up to 5000 * 5000), I am hoping that there is a tool or package
> in MATHEMATICA to solve this?
> 
> Can you please give some link or reference for this?Thank you very much for
> your inforamtion.
> 

What you want is an eigenvector for the eigenvalue of 1. For such 
matrices this is guaranteed to be the largest eigenvalue. So iterations 
of the power method will generally get you such a vector efficiently.

Below I show an example. First we need to generate Markov matrices.

randStochasticMat[dim_] :=
   Transpose[(# / Map[Apply[Plus,#]&, #])& [
     Table[Random[], {dim}, {dim}]]]

For the steady-state vector iterations we will allow settings for 
maximum number of iterations and number of digits of agreement in the 
last two iterations as termination criteria.

ssVector[mmat_, maxiters_:100, tol_:10] := Module[
   {n, vec, oldvec},
   n = Length[mmat];
   vec = Table[Random[], {n}];
   vec = vec / Norm[vec];
   NestWhile [With[{nvec=mat.#}, nvec/Norm[nvec]]&,
     vec, (Max[Abs[#2-#1]]>10^(-tol))&, 2, maxiters]
   ]

Here is an example (run on a 2.8 GHz machine with a Mathematica 
development kernel).

In[4]:= Timing[mat = randStochasticMat[2000];]
Out[4]= {1.7 Second, Null}

In[7]:= Timing[ssvec = ssVector[mat];]
Out[7]= {0.05 Second, Null}

In[10]:= InputForm[Max[Abs[mat.ssvec - ssvec]]]
Out[10]//InputForm= 2.3422583317334045*^-13


Daniel Lichtblau
Wolfram Research



  • Prev by Date: Re: How to include images and graphics in Mathematica
  • Next by Date: Re: Garbage collection problem
  • Previous by thread: any tool/package to solve markov chain?
  • Next by thread: Corner solutions