I have a discrete- time and space markov process that I want to evolve given an initial distribution.

To that, I defined

TimeP[initial_, state_, 0] := initial[state];
TimeP[initial_, state_, t_] := TimeP[initial, state, t] = TimeP[initial, state, t - 1] + ...

where the initial distribution would be something like this:

initial[state_] := 0
initial[firstState] := 1

and ... are probabilities to jump from a state to another state.

Now, with the above definition of TimeP, I can correctly evolve the probabilities. However, when I define another initial distribution, say

initial[state_] := 0
initial[firstState1] := 1

The function TimeP gets confused because it only stores the name initial, and not its definition (i.e. even if initial has now a different definition, TimeP returns the computed value from initial's definition of firstState, and not firstState1.

So, my question is: is a way of telling TimeP that its argument's name doesn't matter, only the actual definition that matters? Or, do you know a better way of implementing this?

This is what I would expect: the function initial is defining a distribution, and the distribution's definition is on its probabilities. When I change initial's definition to another values, TimeP would give another result.

(I'm not using Mathematica's DiscreteMarkovProcess because the transition matrix is almost diagonal with most entries equal to 0, and the number of states is large).

  • 1
    $\begingroup$ It would seem strange to me to rename your initial state whenever you want to specify a new initial distribution (firstState vs firstState1). Can you not simply assign a different value to initial[firstState]? Also, do you have a value set for firstState? $\endgroup$ Oct 21, 2013 at 10:48

1 Answer 1


It is not a direct answer to your question, but it can helps you.

I think it is more efficient to work with sparse arrays and store the resulting distribution in an array, not in DownValues.

n = 100;

m = SparseArray[{Band[{1, 1}] -> 0.5, Band[{1, 2}] -> 0.25, 
   Band[{2, 1}] -> 0.25}, {n, n}];

init = 50;

p0 = MapAt[1. &, #, init] &@ConstantArray[0., n];

p = NestList[m.# &, p0, 100];

ListPlot3D[p, PlotRange -> {0, 0.2}]

enter image description here

Approach with the memoization:

distr[init_, t_] := distr[init, t] = m.distr[init, t - 1];
distr[init_, 0] := m.MapAt[1. &, #, init] &@ConstantArray[0., n];

ListLinePlot[{distr[50, 100], distr[75, 100], distr[50, 200]}]

enter image description here


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.