# How can I simulate the stationary distribution of particles that each moves differently?

Suppose a particle enters a system at $$0.5$$ in the unit interval $$[0,1]$$.

With some probability $$\lambda_{right}$$, particles go right by

$$x_{right} = \frac{x\pi_{H}}{x\pi_{H} + (1-x)\pi_{L} }$$ and with some probability $$\lambda_{left}$$, they go left by $$x_{left} = \frac{x(1-\pi_{H})}{x(1-\pi_{H}) + (1-x)(1-\pi_{L}) }$$

where $$1>\pi_{H}>\pi_{L}>0$$, so that $$x_{right} \geq x \geq x_{left}$$. These are Bayesian updating formulae.

I want to simulate the stationary distribution of this system, where the fraction of particles at each point does not change anymore. At the moment, I don't want to impose any restriction on $$\pi_{H}$$ and $$\pi_{L}$$. If $$\pi_{H} = 1-\pi_{L}$$, for instance, I can simplify the position of each particle just by how many net right moves it had and get closed form solution from a second-order recurrence equation, but this is not I want to do.

(1) Make regular grids from 0 to 1. Depending on $$\pi_{H}$$ and $$\pi_{L}$$, the grid points might not equal the support of positions created in this system, meaning some grid points might not be reached just because of the parameters in the formulae. I just linearly interpolate while finding a fixed point of $$v(x) = \lambda_{right} v(x')+ \lambda_{left} v(x'')$$ such that $$x'_{right} =x$$ and $$x''_{left} = x$$. However, I'm not sure if this is a mathematically or numerically rigorous method. Most importantly, when I impose $$\pi_{H} =1-\pi_{L}$$, it doesn't give me the same simulation result as the closed-form solution.

(2) I think this is a more brut-force way. I make every combination of $$(n,m)$$ where each represents the number of right and left move. The problem is, in this case, I have no clue to what extend I should allow the two natural numbers to be.

Any suggestion or reference would be greatly helpful.

• I think in general, you have to consider all combinations of $(n,m)$ unless the values of $\pi_H$ and $\pi_L$ are such that the system repeats itself. It can be shown that the system has a posibility to repeat a previous configuration iff there is a solution $(k_1,k_2)$ in natural numbers for the equation $(\frac{\pi_H}{\pi_L})^{k_1} = (\frac{1-\pi_L}{1-\pi_H})^{k_2}$ in which case you can remove all (n,m) which are $> (k_1,k_2)$ Feb 23 at 12:25
• Do you mean right/left by or right/left to, because that makes a big difference.
– Ben
Mar 3 at 0:49
• The system is defined on the unit interval, so the updating formulas work without truncating only if I read "by" as "according to" so that $x_{right}$ or $x_{left}$ is the new position. Mar 3 at 11:34

The distribution of $$x$$ can be approximated with a logit-normal distribution for a large number of steps. The distribution will concentrate at 0 or 1, depending on $$\mu$$:

$$\mu=(1-\lambda)\ln\bigg(\frac{1-\pi_H}{1-\pi_L}\bigg)+\lambda\ln\bigg(\frac{\pi_H}{\pi_L}\bigg)$$

• $$\mu<0$$: $$x\rightarrow1$$
• $$\mu>0$$: $$x\rightarrow0$$
• $$\mu=0$$: $$P(x\rightarrow0)=0.5$$; $$P(x\rightarrow1)=0.5$$

I'll assume $$\lambda_{right}+\lambda_{left}=1$$ and use $$\lambda\equiv\lambda_{right}$$.

Let $$x_{a,b}$$ denote the position of a particle that has moved a total of $$a$$ times to the left and $$b$$ times to the right. First, notice that the sequence of the $$n=a+b$$ moves does not affect the value of $$x_{a,b}$$ (e.g., left-left-right-right results in the same position as left-right-right-left).

$$x_{a,b}=\bigg[1+\bigg(\frac{1-\pi_L}{1-\pi_H}\bigg)^a\bigg(\frac{\pi_L}{\pi_H}\bigg)^b\bigg]^{-1}$$

A logit transformation on $$x_{a,b}$$ results in:

$$y_{a,b}=a\ln\bigg(\frac{1-\pi_H}{1-\pi_L}\bigg)+b\ln\bigg(\frac{\pi_H}{\pi_L}\bigg)$$

Possible values of $$y_{a,b}$$ are equally spaced from $$n\ln\Big(\frac{1-\pi_H}{1-\pi_L}\Big)$$ to $$n\ln\Big(\frac{\pi_H}{\pi_L}\Big)$$ with intervals:

$$d=\ln\Big(\frac{\pi_H}{\pi_L}\Big)-\ln\Big(\frac{1-\pi_H}{1-\pi_L}\Big)$$

$$b\sim\text{Bin}(n,\lambda)$$, so $$y_{a,b}$$ follows a scaled and shifted binomial distribution:

$$b=\frac{y_{a,b}-n\ln\Big(\frac{1-\pi_H}{1-\pi_L}\Big)}{n\cdot{d}}\sim\text{Bin}(n,\lambda)$$

We can use the normal approximation to the binomial as $$n$$ gets large:

$$y_{a,b}\sim\mathcal{N}(n\mu,d^2n\lambda(1-\lambda))$$

Equivalently, $$x_{a,b}$$ can be approximated with the logit-normal distribution.

As $$n$$ increases, the magnitude of the ratio of the mean to the standard deviation increases if $$\mu\neq0$$, so $$x$$ will tend to either 0 or 1. For $$\mu=0$$, the mean of $$x$$ remains at 0.5, and, although the most probable positions are near 0.5, they make up an increasingly smaller proportion of the total probability.