Re: Distribution of state occupancies in a multistate Markov model

*To*: mathgroup at smc.vnet.net*Subject*: [mg104449] Re: Distribution of state occupancies in a multistate Markov model*From*: mathlawguy <chandler.seth at gmail.com>*Date*: Sun, 1 Nov 2009 03:57:27 -0500 (EST)*References*: <hc92bk$edu$1@smc.vnet.net> <hce3vp$r1u$1@smc.vnet.net>

On Oct 30, 2:18 am, fd <fdi... at gmail.com> wrote: > Seth. > > Well, I'm not sure if my approach is the one that suits what you need, > but let me try. > > I would start looking if the chapman-kolmogorov equation can be > analytically solved > > P(x1,t1|x0,t0)=Integrate P(x1,t1|x2,t2)P(x2,t2|x0,t0) dx2 > > integrate over all possible states x2. Depending on your transition > matrix, which are the P(x1,t1|x2,t2)'s, you may be able to perform > this integration.. > > please, take a look in Gardiner's book I mention below..on mackay's > book there are a few comments on chap 30.. > > refs: > > David J. C. MacKay, Information Theory, Inference and learning > algorithms, Cambridge > > You can find a copy of this book on-line if you agree not to print > it.. > > another good reference is the book by the new zealander physicist, C. > W. Gardiner, Handbook of Stochastic methods > > good luck Thanks for the ideas. I will definitely check out the references you've provided. In the meantime, I think I have made some definite headway on my problem. Let P be the sequence of possibly heterogeneous stochastic matrices. We can define Phi as the sequence of their dot products. \[CapitalPhi][P:{___?MatrixQ}]:=FoldList[Dot,IdentityMatrix[Dimensions [P[[1]]][[1]]],P] I believe Phi[[n]][[i,j]] now represents the probability that if you start at state i you will be in state j at time n OK. So we now define a sequence of Bernoulli random variables with probability Phi[[m]][[i,j]] as m runs from 1 to n. If we take the convolution of these Bernoulli random variables, we get the probability distribution for the number of times we will have been in state j having started at state i and having made a total of n transitions. In Mathematica, it appears that one can do this by taking the CharacteristicFunction of the Bernoulli random variables, multiplying them and then taking the Inverse Fourier Transform with FourierParameters->{1,1}. Here's the actual code I'm using. dist[P : {___?MatrixQ}, {i_Integer, j_Integer}]:=InverseFourierTransform[ Times @@ Map[ CharacteristicFunction[BernoulliDistribution[#[[i, j]]], t] &, \[CapitalPhi][P]] // Expand, t, x, FourierParameters -> {1, 1}] When I do this, I end up with an expression involving a bunch of DiracDeltas, which is what one would expect. Moreover the process appears pretty fast. I've develop an analytic expression for the coefficients in front of the DiracDeltas involving binomials but for some reason it takes longer to compute them than it does to just use the method described above. Moreover, if one wants to compute the probability distribution of the number of times the process will have been in some set of states J, one sets the parameters of the Bernoulli random variables as Sum[Phi [[n]][i,j]], with j an element of J and then performs a similar computation. Here's the code: dist[P : {___?MatrixQ}, {i_Integer, j : {__Integer}}] := InverseFourierTransform[ Times @@ Map[ CharacteristicFunction[BernoulliDistribution[Total[#[[i, j]]]], t] &, \[CapitalPhi][P]] // Expand, t, x, FourierParameters -> {1, 1}] If anyone thinks I am doing this incorrectly or has a better/faster idea, please feel strongly encouraged to let me know as I am about to base some important computations on the idea.