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

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.
