In a model I have a discrete two-state first order Markov process, defined by a (2x2) transition matrix with two free parameters. If the first state occurs then the process outputs zero for that period. If the second state occurs then the process outputs a Gamma-distributed random variable characterised by another two parameters.

I can work out the theoretical mean and variance of the overall stationary process using pen and paper, but I am unsure how to check my results (and also I want to practice my Mathematica) using symbolic calculations. Put another way, I would like to know how to calculate the mean and variance to which the composite Markov-Gamma process will converge. My code is below - all I have currently is a way to calculate the mean of the Markov process, P, and its stationary mean (Mean[S]). However, I would like to be able to derive the mean and variance of R inserted into P for times when the process is in state 2, and zeros for when P is in state 1. Many thanks! Ben

M = {{1 - pdw, pdw}, {pwd, 1 - pwd}}
P  = DiscreteMarkovProcess[{1, 0}, M]
S = StationaryDistribution[P]
R = GammaDistribution[a,b]
  • $\begingroup$ things like PDF[P[n], k] // PiecewiseExpand or Mean[P[n]] // FullSimplify? $\endgroup$
    – chris
    Aug 21, 2014 at 17:09
  • $\begingroup$ I did not understand the link with R? $\endgroup$
    – chris
    Aug 21, 2014 at 17:23
  • $\begingroup$ R is a gamma distribution which determines the outcome of the series if it is in state '2'. Otherwise the process outputs zero in state '1'. $\endgroup$
    – ben18785
    Aug 21, 2014 at 17:46

2 Answers 2



$Assumptions = {0 <= pdw <= 1, 0 <= pwd <= 1};
trM = {{1 - pdw, pdw}, {pwd, 1 - pwd}};
dmP = DiscreteMarkovProcess[{1, 0}, trM];
gD = GammaDistribution[a, b];

the mean and variance of the described process for a time slice at time t are

mean = Expectation[(m[t] - 1)*g, {m \[Distributed] dmP, g \[Distributed] gD}] 
variance = Expectation[((m[t] - 1)*g - mean)^2, {m \[Distributed] dmP, g \[Distributed] gD}] // Simplify
-((a b pdw (-1 + (1 - pdw - pwd)^t))/(pdw + pwd))
-((a b^2 pdw (-1 + (1 - pdw - pwd)^t) (pdw + a pdw (1 - pdw - pwd)^t + pwd + a pwd))/(pdw + pwd)^2)

The stationary mean and variance are

stationaryMean = Limit[mean, t -> ∞]
stationaryVariance = Limit[variance, t -> ∞]
(a b pdw)/(pdw + pwd)
(a b^2 pdw (pdw + pwd + a pwd))/(pdw + pwd)^2

An alternative approach is to define a TransformedDistribution

trD = TransformedDistribution[(m - 1)* g, {m \[Distributed] SliceDistribution[dmP, t], 
       g \[Distributed] gD}]

and then calculate the stationary mean and variance with

Mean@trD ~Limit~ (t -> ∞)
Variance@trD ~Limit~ (t -> ∞)

Performing simulations:
One can use the TransformedDistribution to perform simulations:

rvSim = RandomVariate[trD /. {pdw -> 0.2, pwd -> 0.7, a -> 1, b -> 2, t -> #}] & /@ Range[1000];

Nevertheless, the calculation time is much lower, if one first uses RandomFunction with the DiscreteMarkovProcess and than does the state dependent transformation:

rf = RandomFunction[dmP /. {pdw -> 0.2, pwd -> 0.7}, {1, 10000}];

f[x_Integer] := 0 /; x == 1
f[x_Integer] := RandomVariate[GammaDistribution[1, 2]] /; x == 2

sim = f /@ (rf["Path"][[All, 2]]);
  • $\begingroup$ I am impressed you seem to have understood what the OP wanted based on limited evidence :-) $\endgroup$
    – chris
    Aug 22, 2014 at 7:44
  • $\begingroup$ Thanks to Karsten for answering my question - that helps me understand how to approach this problem. $\endgroup$
    – ben18785
    Aug 22, 2014 at 9:23
  • $\begingroup$ @Karsten. Can I ask how I can simulate draws from the transformed distribution? I have tried RandomFunction (substituting in the required parameters with numbers), but this does not seem to work. $\endgroup$
    – ben18785
    Aug 22, 2014 at 10:03
  • $\begingroup$ @user2003951 your difficulties seem to be due to the circumstance that for Mathematica there is a bigger difference between distributions and random processes than for you. RandomFunction is for processes, not for distributions. For distributions you should use RandomVariate. I'll make an edit to my answer, in order to give you an example. $\endgroup$
    – Karsten7
    Aug 22, 2014 at 11:09
  • $\begingroup$ @Karsten - that's great thank you. Yes, the latter option is much faster. Can I ask though, how do I actually plot the data from the simulation? I have tried ListPlot[sim] but that doesn't seem to work. Sorry! $\endgroup$
    – ben18785
    Aug 22, 2014 at 13:13

Following closely the documentation for DiscreteMarkovProcess (I am no expert though…)

M = {{1 - w, w}, {d, 1 - d}} P = DiscreteMarkovProcess[{1, 0}, M]

Mathematica graphics

Mean[P[n]] // FullSimplify

Mathematica graphics

Variance[P[n]] // FullSimplify

Mathematica graphics


Mathematica graphics

data = RandomFunction[P /. w -> 1/2 /. d -> 1/4, {0, 10}]

Mathematica graphics

ListPlot[data, Filling -> Axis, Ticks -> {Automatic, {1, 2, 3}}]

Mathematica graphics

MarkovProcessProperties[P /. w -> 1/2 /. d -> 1/4]

Mathematica graphics

Now for the stationary process

S = StationaryDistribution[P] // FullSimplify;

PDF[S, n]

Mathematica graphics

Mean[S] // FullSimplify

Mathematica graphics

Variance[S] // FullSimplify

Mathematica graphics

  • $\begingroup$ Hi, Thanks for your answer, although it is not quite what I am looking for. I can derive the properties of a simple Markov chain using Mathematica. What I can't do is form the composite process consisting of zeros if the process is in state 1, and a gamma distributed random variable in state 2. I don't know how to make this composite process using Mathematica, then symbollically derive its long run mean and variance. Thanks! $\endgroup$
    – ben18785
    Aug 21, 2014 at 17:50
  • $\begingroup$ ok I cannot help you there because I know too little about Markov processes. There are ways to compose probabilities in mathematica: may be someone else will help. $\endgroup$
    – chris
    Aug 21, 2014 at 18:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.