Let $\mu_0 (dx)$ be the standard one-dimensional Cauchy distribution, i.e.
\begin{align} \mu_0 (dx) = \frac{1}{\pi} \frac{1}{1+x^2} dx. \end{align}
Suppose I fix $h \in [0, 1]$, and form a Markov chain $\{X_n\}_{n \geqslant 0}$ as follows:
- At step $n$, I sample $Y_n \sim \mu_0$.
- I then set $X_{n+1} = (1 - h) X_n + h Y_n$
It is not so hard to show that this chain admits $\mu_0$ as a stationary measure, as this essentially comes from the fact that Cauchy distribution is a stable distribution.
What I'm interested in is the reversal of this Markov chain. More precisely, if the chain I describe above uses the Markov kernel $q (x \to dy)$, I want to understand the Markov kernel $r$ such that
\begin{align} \mu_0 (dx) q (x \to dy) = \mu_0 (dy) r (y \to dx). \end{align}
Fortunately, all of the quantities involved have densities with respect to Lebesgue measure, and as such, I can write down what $r (y \to dx)$ is:
\begin{align} r (y \to dx) = \frac{1 + y^2}{\pi} \frac{1}{1 + x^2} \frac{h}{ h^2 + (y - (1 - h) x)^2} dx. \end{align}
My question is then: is there a simple, elegant way to draw exact samples from $r$?
I would highlight that this is not a purely algorithmic question; I'd really like to understand what this reversal kernel $r$ is doing. A nice byproduct of that would then be that I could simulate from it easily.
For completeness, some of the `purely algorithmic' solutions I had considered were the following.
- I could try rejection sampling, and in principle this would work, but it wouldn't really give me insight into the nature of the Markov chain.
- I could try something like the inverse CDF method, but it seems to me that the CDF of $r$ is not particularly nice to work with. As such, I'd have to use e.g. Newton iterations to use this method, and I'd prefer to not have to do this.
