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

Question

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.

I had two options.

(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.

Answer 1

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.