# Calculating probability of results of a non-deterministic function

I am working on a computer program. It can be represented as a linear series of states $s_0 \ldots s_n$

Function $N(s_x)$ is defined such that:

• With a probability $p$ it will return the next state $$s_x \rightarrow \left\{ \begin{array}{ll} s_{x+1} & \mbox{if } x \lt n \\ s_{x} & \mbox{otherwise} \end{array} \right.$$
• With a probability of $p$ it will return the previous state $$s_x \rightarrow \left\{ \begin{array}{ll} s_{x-1} & \mbox{if } x \gt 0 \\ s_{x} & \mbox{otherwise} \end{array} \right.$$
• With a probability of $1-2p$ it will return the current state $$s_x \rightarrow s_x$$

The program is an iterative application of $N$ to the initial value $s_i$ $r$ times $$Program = N^r(s_i)$$

I need to find $P(N^r(s_i) = s_x)$ for all $0 \leq x \leq n$ (the probability that applying $N$ iteratively to $s_i$ $n$ times yields the state $s_x$)

For small values I've been manually working out the different combinations that can lead to a particular state. For example, if $n=2$, $i=1$, $p=0.25$ and $r=2$ I can create a tree to help figure out the probabilities: In this case: $$P(s_0) = 0.25*0.25 + 0.25*0.5 + 0.5*0.25 = 0.3125$$ $$P(s_1) = 0.25*0.25 + 0.5*0.5 + 0.25*0.25 = 0.375$$ $$P(s_2) = 0.5*0.25 + 0.25*0.5 + 0.25*0.25 = 0.3125$$

However, this obviously gets impractical for the larger values I need to deal with in real life (e.g. $r=500$). Is there an easier way to calculate the probabilities?

## 3 Answers

You've got a time-homogeneous Markov chain with finite state space, and you can diagonalize the transition matrix and decompose your initial state into its eigenvectors; then applying $N$ $r$ times just becomes multiplying each component with the $r$-th power of the corresponding eigenvalue. Of course this will only be more efficient than computing the $r$ steps directly if $n \ll r$.

Your probability is flowing according to the heat equation, and your boundary conditions are equivalent to a mirror, so if $i$=20 and $n$=99, then $s_x$ is the same as for an infinite line of heat flow initialized with point hotspots at ...,-221,-180,-21,20,179,220,279,... .

The heat at a point is just the sum of the contributions from these spots, and the contribution from a spot can be quickly estimated as $\frac{\exp{(\frac{-d^2}{4pr})}}{\sqrt{4\pi pr}}$, where $d$ is the distance from the heat source to the point of measurement. These contributions die off very quickly (although any spot can contribute in principle if the time $r$ gets big enough), so you don't need to include spots from far away in your estimate.

1. Exact Values

You get exact values using matrix multiplication, as explained in joriki's answer: $P(N^r(s_i) = s_x)$ is the value in the $x$th column and $i$th row of the $r$th power of the transition matrix $P$. In your worked out example, the matrix is $$P=\pmatrix{3/4&1/4&0\cr 1/4&1/2&1/4\cr 0&1/4&3/4}.$$ The second power is $$P^2=\pmatrix{10/16&5/16&1/16\cr 5/16&6/16&5/16\cr 1/16&5/16&10/16}$$ and (as expected) your worked out probabilities $.3125$, $.3750$, and $.3125$ are the second row of $P^2$. That's because the second row corresponds to the starting position $i=1$. (Note: The rows and columns of $P$ are labelled from $0$ to $n$.)

For larger $r$, you take higher powers of the matrix. When $r=10$ we get $$P^{10}=\pmatrix{ {189525\over 524288}&{349525\over 1048576}&{320001\over 1048576}\cr {349525\over 1048576}& {174763\over 524288}&{349525\over 1048576}\cr {320001\over 1048576}& {349525\over 1048576}& {189525\over 524288}}$$

2. Approximation

For large $r$ values, $P(N^r(s_i) = s_x)\approx 1/(n+1)$ for any $i$ and $x$. That is because the equilibrium distribution for your process is uniform on the state space $\{0,1,\dots,n\}$. In the long run, your process is equally likely to be in any state, regardless of starting position.

For instance, in your worked out example, the probabilities are all approximately $1/3$. When $r=2$, the approximation isn't that great, but the convergence is exponentially fast. For $r=10$, the middle row of the matrix $P^{10}$ says that $P(s_0)=.3333330154$, $P(s_1)=.3333339691$, and $P(s_2)=.3333330154$. These are all pretty close to $1/3$.