Defining a function of a distribution

Question

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

Clear[initial]
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).

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}]

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]}]