<think> We are given a Markov chain defined by: X_{n+1} = (1-h) * X_n + h * Y_n, where Y_n ~ Cauchy(0,1) (standard Cauchy)

We know that the stationary distribution is the standard Cauchy, i.e., μ_0(dx) = (1/(π(1+x^2))) dx.

The forward kernel q(x -> dy) is the density of the distribution of X_{n+1} given X_n = x. Since: X_{n+1} = (1-h)x + hY, and Y ~ Cauchy(0,1), then the distribution of X_{n+1} given x is Cauchy with location (1-h)x and scale h. Why? Because the Cauchy distribution is stable: scaling and shifting a Cauchy random variable by constants gives another Cauchy. Specifically, if Y ~ Cauchy(0,1), then aY + b ~ Cauchy(b, a). Here, we have: h * Y ~ Cauchy(0, h) [since scaling by h gives scale h] then adding (1-h)*x: (1-h)x + hY ~ Cauchy((1-h)*x, h)

Therefore, the forward kernel density is: q(x -> y) = (1/(π h)) * 1 / [1 + ((y - (1-h)*x)/h)^2] = (1/(π h)) * h^2 / [h^2 + (y - (1-h)*x)^2] = h / [π (h^2 + (y - (1-h)*x)^2)]

Now, the reversal kernel r(y -> dx) must satisfy the detailed balance condition with respect to the stationary measure μ_0: μ_0(dx) q(x -> dy) = μ_0(dy) r(y -> dx)

Since we have densities, we can write: [μ_0(x) dx] * [q(x -> y) dy] = [μ_0(y) dy] * [r(y -> x) dx]

Therefore, solving for r(y -> x): r(y -> x) = [μ_0(x) * q(x -> y)] / μ_0(y)

Plugging in the densities: μ_0(x) = 1/(π(1+x^2)) μ_0(y) = 1/(π(1+y^2)) q(x->y) = h / [π (h^2 + (y - (1-h)*x)^2)]

So: r(y -> x) = [ (1/(π(1+x^2))) * (h/(π (h^2 + (y - (1-h)*x)^2))) ] / [1/(π(1+y^2))] = [ (1/(π(1+x^2))) * (h/(π (h^2 + (y - (1-h)*x)^2))) ] * [π(1+y^2)] = [h (1+y^2)] / [ π (1+x^2) (h^2 + (y - (1-h)*x)^2) ]

This matches the expression given in the question.

Now, the goal is to sample from the distribution with density: r(y -> x) = [h (1+y^2)] / [ π (1+x^2) (h^2 + (y - (1-h)*x)^2) ]

We want to understand this distribution and find an elegant way to sample from it.

Let me rewrite the expression for the density. Notice that the denominator has two terms: (1+x^2) and (h^2 + (y - (1-h)*x)^2). We can interpret this as a product of two Cauchy-like densities.

Consider the term: h^2 + (y - (1-h)*x)^2 = h^2 + (y - (1-h)x)^2

We can rearrange the entire density:

r(y->x) = [h (1+y^2)] / [ π (1+x^2) (h^2 + (y - (1-h)x)^2) ]

Notice that the stationary measure μ_0(x) appears in the denominator (via 1/(1+x^2)), and the forward kernel q(y->x) would be similar to the forward kernel but with the roles of x and y swapped? Not exactly.

Alternatively, let us try to relate this to the forward kernel. The forward kernel from x to y is: q(x->y) = h / [π (h^2 + (y - (1-h)x)^2)]

But the reversal kernel is: r(y->x) = [μ_0(x)/μ_0(y)] * q(x->y)

This is the general form. Since we know that the chain is reversible if and only if the detailed balance holds, and in this case, is it reversible? Only if r(y->x) = q(y->x). But let's check:

The forward kernel from y to x would be: q(y->x) = h / [π (h^2 + (x - (1-h)y)^2)]

But note that: r(y->x) = [μ_0(x)/μ_0(y)] * q(x->y) = [ (1+x^2)^{-1} / (1+y^2)^{-1} ] * q(x->y) = [(1+y^2)/(1+x^2)] * q(x->y)

and q(x->y) = h / [π (h^2 + (y - (1-h)x)^2)]

So unless (1+y^2)/(1+x^2) * q(x->y) = q(y->x), which is not generally true, the chain is not reversible.

But we can try to interpret r(y->x) as a mixture or as a conditional distribution.

Another idea: let us consider the joint distribution of (X_n, X_{n+1}) under stationarity. Since the chain is stationary, the joint density is: p(x,y) = μ_0(x) * q(x->y) = [1/(π(1+x^2))] * [h/(π (h^2 + (y - (1-h)x)^2))]

The reversal kernel r(y->x) is the conditional density of X_n given X_{n+1}=y, because in the reversed chain, we have: P(X_n = x | X_{n+1} = y) = r(y->x)

But note that the joint density of (X_n, X_{n+1}) is the same as the joint density of (X_{n+1}, X_n) for the reversed chain. Therefore, the conditional distribution of X_n given X_{n+1}=y is exactly r(y->x).

So we can write: p(x|y) = p(x,y) / p(y) = [μ_0(x) q(x->y)] / μ_0(y) = r(y->x)

So we are trying to sample from the conditional distribution of X_n given X_{n+1}=y, under stationarity.

Now, how can we sample from this conditional distribution?

Recall the model: X_{n+1} = (1-h) * X_n + h * Y_n, with Y_n ~ Cauchy(0,1) independent of X_n.

But under stationarity, X_n ~ Cauchy(0,1). Therefore, we have a linear relationship: X_{n+1} = (1-h) * X_n + h * Y_n.

Given X_{n+1}=y, we have: y = (1-h) * X_n + h * Y_n.

We know that X_n and Y_n are independent, and both are Cauchy(0,1). Therefore, we are trying to sample X_n given that (1-h)X_n + hY_n = y.

This is a linear constraint. Since both X_n and Y_n are Cauchy, we can use properties of the Cauchy distribution.

Note: The sum of two independent Cauchy random variables is Cauchy, but here we have a constraint that they are linearly combined to a fixed value.

We can think of the conditional distribution of X_n given that: (1-h) * X_n + h * Y_n = y.

Since Y_n = [y - (1-h)*X_n] / h, and Y_n must be a standard Cauchy.

Therefore, the conditional density of X_n given X_{n+1}=y is proportional to the joint density of (X_n, Y_n) but with the constraint. The joint density of (X_n, Y_n) is: f_{X_n, Y_n}(x, z) = [1/(π(1+x^2))] * [1/(π(1+z^2))]

But we set z = [y - (1-h)*x] / h.

Therefore, the conditional density of X_n given X_{n+1}=y is proportional to: f_{X_n}(x) * f_{Y_n}( [y - (1-h)*x] / h ) = [1/(π(1+x^2))] * [1/(π(1 + ([y - (1-h)*x]/h)^2))]

But note that we have to multiply by the Jacobian? Actually, in the conditional density, we are not changing variables; we are evaluating the joint density at (x, [y-(1-h)x]/h) and then normalizing. However, the normalizing constant is the marginal density of X_{n+1} at y, which we know is μ_0(y) = 1/(π(1+y^2)). But we already have the expression for r(y->x) which is the conditional density.

So we have: f_{X_n|X_{n+1}=y}(x) ∝ 1/[(1+x^2) * (1 + ( (y - (1-h)x)/h )^2 )]

But note that: 1 + ( (y - (1-h)x)/h )^2 = [h^2 + (y - (1-h)x)^2] / h^2

So: f_{X_n|X_{n+1}=y}(x) ∝ 1/[(1+x^2) * (h^2 + (y - (1-h)x)^2) / h^2] ∝ h^2 / [(1+x^2) (h^2 + (y - (1-h)x)^2)]

But wait, this is proportional to: h^2 / [(1+x^2) (h^2 + (y - (1-h)x)^2)]

However, we know from the detailed balance that the conditional density is: r(y->x) = [h (1+y^2)] / [ π (1+x^2) (h^2 + (y - (1-h)x)^2) ]

So indeed, the conditional density is proportional to 1/[(1+x^2) (h^2 + (y - (1-h)x)^2)], and the constant is [h (1+y^2)]/(π) times the normalizing constant? Actually, we see that the expression we have for r(y->x) is exactly the conditional density.

Now, to sample from this distribution, we can use the representation in terms of the latent variable. Notice that:

y = (1-h) * x + h * z, where z is an independent standard Cauchy.

Therefore, given y, we have: z = (y - (1-h)*x) / h.

But since z must be a standard Cauchy, the conditional density of x given y is as above.

Alternatively, we can use the fact that the joint distribution of (X_n, X_{n+1}) is a bivariate Cauchy? Not exactly, but we can use the representation with latent variables.

Another idea: we can express the conditional distribution as a mixture or as a scale mixture of normals? But Cauchy distributions are stable and have representations as normal scale mixtures.

Recall that a standard Cauchy can be represented as: X = Z / sqrt(W), where Z ~ N(0,1) and W ~ Gamma(1/2, 1/2) (i.e., chi-square with 1 df) independent.

But here we have a more complex structure.

Alternatively, we can use the fact that the ratio of two independent standard normals is Cauchy. But that might not help directly.

Let me try to rewrite the conditional density:

r(y->x) ∝ 1/[(1+x^2) * (h^2 + (y - (1-h)x)^2)]

We can decompose the denominator. Consider the quadratic form in x.

Expand the denominator:

(1+x^2) * [h^2 + (y - (1-h)x)^2] = (1+x^2) * [h^2 + y^2 - 2y(1-h)x + (1-h)^2 x^2] = (1+x^2) * [ (1-h)^2 x^2 - 2y(1-h)x + (y^2 + h^2) ]

But this is a quartic in x, which is messy.

Alternatively, we can complete the square in the exponent? But the density is rational, not exponential.

Another approach: use the fact that the Cauchy distribution is closed under linear fractional transformations. In fact, the conditional distribution might be a Cauchy distribution or a mixture of Cauchys.

Let me assume that the conditional distribution is a Cauchy distribution. Then it would be characterized by its location and scale. But the expression for the density is a ratio of a quadratic in the numerator (which is constant in x) and a quartic in the denominator? Actually, the numerator of r(y->x) is constant in x? No, the numerator has (1+y^2) which is constant for fixed y, but the denominator is (1+x^2) times a quadratic in x.

So the density is proportional to: 1 / [ (1+x^2) * (A x^2 + B x + C) ]

where A = (1-h)^2, B = -2y(1-h), C = h^2 + y^2.

This is a rational function that can be decomposed by partial fractions. Let me try that.

We want to write: 1 / [ (1+x^2) * (A x^2 + B x + C) ] = [D x + E]/(1+x^2) + [F x + G]/(A x^2 + B x + C)

But note that the denominator factors are quadratic, so we can write:

1 / [ (1+x^2) (A x^2 + B x + C) ] = [P x + Q]/(1+x^2) + [R x + S]/(A x^2 + B x + C)

Multiplying both sides by (1+x^2)(A x^2 + B x + C) and equating coefficients, we get a system of linear equations.

However, this might be messy, but it is possible. Then we would have the density as a sum of two terms, each of which is proportional to a Cauchy-like density?

But note: the term [P x + Q]/(1+x^2) is proportional to the density of a Cauchy distribution if we can write it as a constant times 1/(1+x^2) plus a constant times x/(1+x^2). However, the term x/(1+x^2) is the derivative of the log of the Cauchy density? Actually, it is the negative of the derivative of the Cauchy density. But that doesn't give a probability density.

Alternatively, we can recognize that the entire expression might be represented as a mixture of two Cauchy distributions.

Let me consider the partial fraction decomposition. Set:

[P x + Q]/(1+x^2) + [R x + S]/(A x^2 + B x + C) = 1 / [ (1+x^2) (A x^2 + B x + C) ]

Then:

(P x + Q)(A x^2 + B x + C) + (R x + S)(1+x^2) = 1

Expand:

P x (A x^2 + B x + C) + Q (A x^2 + B x + C) + R x (1+x^2) + S (1+x^2) = 1

= P (A x^3 + B x^2 + C x) + Q (A x^2 + B x + C) + R (x + x^3) + S (1 + x^2)

= (A P + R) x^3 + (B P + A Q + S) x^2 + (C P + B Q + R) x + (C Q + S)

Set equal to 1, which is 0x^3 + 0x^2 + 0*x + 1.

So we have the system:

A P + R = 0 B P + A Q + S = 0 C P + B Q + R = 0 C Q + S = 1

This is a linear system in P, Q, R, S.

But note: A = (1-h)^2, B = -2y(1-h), C = h^2 + y^2.

This system is solvable, but the expressions might be complicated. However, we can solve it.

From the first equation: R = -A P.

From the fourth equation: S = 1 - C Q.

Plug into the second and third equations.

Second equation: B P + A Q + (1 - C Q) = 0 => B P + A Q - C Q = -1 => B P + (A - C) Q = -1.

Third equation: C P + B Q + R = C P + B Q - A P = (C - A) P + B Q = 0.

So we have:

(C - A) P + B Q = 0 ... (eq3) B P + (A - C) Q = -1 ... (eq2)

But note that A - C = - (C - A), so eq2 becomes: B P - (C - A) Q = -1.

Set D = C - A, then we have:

D P + B Q = 0 ... (eq3) B P - D Q = -1 ... (eq2')

Now, we can write:

From eq3: P = - (B / D) Q.

Plug into eq2': B [ - (B / D) Q ] - D Q = -1 => - (B^2 / D) Q - D Q = -1 => Q [ - (B^2 + D^2) / D ] = -1 => Q (B^2 + D^2) / D = 1 => Q = D / (B^2 + D^2)

Then P = - (B / D) * [D / (B^2 + D^2)] = -B / (B^2 + D^2)

Now, D = C - A, so:

P = -B / (B^2 + (C-A)^2) Q = (C-A) / (B^2 + (C-A)^2)

Then R = -A P = A B / (B^2 + (C-A)^2) S = 1 - C Q = 1 - C (C-A) / (B^2 + (C-A)^2)

But we don't really need R and S for the decomposition of the density? Actually, we only need the first part for the Cauchy component? Not exactly, because the density is then:

r(y->x) ∝ [P x + Q]/(1+x^2) + [R x + S]/(A x^2 + B x + C)

But this is not necessarily a mixture of two densities. It is a sum, and it might be negative in parts? We need to check the signs.

Note that the denominator is always positive, so the entire expression is positive. But the numerators P x + Q and R x + S might be negative for some x. However, when we integrate, the total integral is positive. But for sampling, we want a representation as a mixture of two distributions that we can sample from.

Alternatively, we can note that the expression:

1 / [ (1+x^2) (A x^2 + B x + C) ] = [P x + Q]/(1+x^2) + [R x + S]/(A x^2 + B x + C)

and then the density r(y->x) is proportional to this. Therefore, we can write:

r(y->x) = K(y) * { [P x + Q]/(1+x^2) + [R x + S]/(A x^2 + B x + C) }

where K(y) is a normalizing constant that depends on y.

But note that the term [P x + Q]/(1+x^2) is not a probability density by itself because P x + Q might be negative. Similarly for the other term.

However, we can try to write it as a mixture if we can show that both parts are non-negative and integrate to positive numbers. But that might not be the case.

Another idea: we can use the representation with latent variables to sample from the conditional distribution.

Recall that we have: y = (1-h) * x + h * z, where x and z are independent standard Cauchy.

Given y, we want to sample x such that the equation holds for some z that is standard Cauchy.

We can solve for z: z = (y - (1-h)*x) / h.

The constraint is that z must be a standard Cauchy. Therefore, the conditional density of x given y is proportional to the product of the density of x and the density of z evaluated at (y - (1-h)*x)/h.

That is: f(x|y) ∝ f_X(x) * f_Z( (y - (1-h)*x)/h ) = [1/(π(1+x^2))] * [1/(π(1 + [(y - (1-h)*x)/h]^2))]

which is exactly what we have.

Now, we can use a latent variable approach to sample from this distribution. Specifically, we can use a Gibbs sampler or a slice sampler, but that would be approximate and not exact.

Alternatively, we can recognize that the conditional distribution is a Cauchy distribution with parameters that depend on y.

Let me assume that the conditional distribution of x given y is Cauchy(location=m, scale=s) for some m and s that depend on y. Then the density would be:

f(x|y) = 1 / [π s (1 + ((x-m)/s)^2)] = s / [π (s^2 + (x-m)^2)]

But our density is:

r(y->x) = [h (1+y^2)] / [ π (1+x^2) (h^2 + (y - (1-h)x)^2) ]

We can try to match the two expressions. However, the denominator has two factors: (1+x^2) and (h^2 + (y - (1-h)x)^2), while the Cauchy density has only one quadratic factor. So it is not a Cauchy distribution.

But perhaps it is a mixture of two Cauchy distributions? Let me go back to the partial fractions.

From earlier, we have:

f(x|y) ∝ 1 / [ (1+x^2) (A x^2 + B x + C) ] = [P x + Q]/(1+x^2) + [R x + S]/(A x^2 + B x + C)

Now, if we can write this as a mixture, then we might be able to sample by choosing one of the two components with appropriate probability and then sampling from a Cauchy or from the other distribution.

But the other distribution is a scaled and shifted Cauchy? Note that the second denominator A x^2 + B x + C is a quadratic that may or may not be written as a constant times (1 + ((x-m)/s)^2). In general, it can be written as:

A x^2 + B x + C = A [ x^2 + (B/A) x + C/A ] = A [ (x + B/(2A))^2 + (C/A - B^2/(4A^2)) ]

So if the discriminant is positive, then it is a constant times (1 + ((x - m)/s)^2) for some m and s. But the discriminant is:

D = B^2 - 4A C = [-2y(1-h)]^2 - 4*(1-h)^2*(h^2+y^2) = 4y^2 (1-h)^2 - 4 (1-h)^2 (h^2+y^2) = 4 (1-h)^2 [ y^2 - (h^2+y^2) ] = 4 (1-h)^2 [ -h^2 ] = -4 (1-h)^2 h^2 < 0

So the discriminant is negative, which means that the quadratic A x^2 + B x + C is always positive and has no real roots, and it is always at least:

A * [ (x + B/(2A))^2 + (4AC - B^2)/(4A^2) ] = A * [ (x + B/(2A))^2 + ( - ( -4A C + B^2) ) / (4A^2) ] = A * [ (x + B/(2A))^2 + (4A C - B^2)/(4A^2) ]

But since 4A C - B^2 = 4*(1-h)^2*(h^2+y^2) - [ -2y(1-h) ]^2 = 4(1-h)^2 [ (h^2+y^2) - y^2 ] = 4(1-h)^2 h^2 > 0.

So we can write:

A x^2 + B x + C = A [ (x + B/(2A))^2 + (4A C - B^2)/(4A^2) ] = A [ (x - m)^2 + s^2 ]

where m = -B/(2A) = [2y(1-h)] / [2 (1-h)^2] = y / (1-h) [since 1-h != 0, and if h=1 then we handle separately]

and s^2 = (4A C - B^2) / (4A^2) = [4 (1-h)^2 h^2] / [4 (1-h)^4] = h^2 / (1-h)^2.

So s = |h / (1-h)|, and since h in [0,1], s = h/(1-h) if h<1, and if h=1, then we have to handle separately.

So for h in [0,1) and h>0, we have:

A x^2 + B x + C = (1-h)^2 [ (x - y/(1-h))^2 + (h/(1-h))^2 ]

Therefore, the second term in the partial fraction decomposition is:

[R x + S] / [ (1-h)^2 ( (x - y/(1-h))^2 + (h/(1-h))^2 ) ]

This is proportional to a Cauchy density with location y/(1-h) and scale h/(1-h), but multiplied by (R x + S).

Specifically, the density of a Cauchy(m, s) is:

g(x) = 1 / [π s (1 + ((x-m)/s)^2)] = s / [π (s^2 + (x-m)^2)]

So here, the denominator is (1-h)^2 [ (x - m)^2 + s^2 ] with m = y/(1-h), s = h/(1-h), so:

s^2 = (h/(1-h))^2.

Therefore, the term:

1 / (A x^2 + B x + C) = 1 / [ (1-h)^2 ( (x - m)^2 + s^2 ) ] = [1/( (1-h)^2 )] * 1/( (x-m)^2 + s^2 )

But the Cauchy density would have s / [π ( (x-m)^2 + s^2 ) ], so:

1/( (x-m)^2 + s^2 ) = [π s]^{-1} * [s / (π ( (x-m)^2 + s^2 ))] * π s = [π s]^{-1} * g(x) * π s

This is not directly helpful.

Now, let's go back to the partial fraction decomposition. We have:

1 / [ (1+x^2) (A x^2 + B x + C) ] = [P x + Q]/(1+x^2) + [R x + S]/(A x^2 + B x + C)

We have expressions for P and Q in terms of B and D = C - A:

P = -B / (B^2 + D^2) Q = D / (B^2 + D^2)

where D = C - A.

Now, A = (1-h)^2, B = -2y(1-h), C = h^2 + y^2.

So D = C - A = (h^2 + y^2) - (1-h)^2 = h^2 + y^2 - (1 - 2h + h^2) = y^2 + 2h - 1

Then B^2 = [-2y(1-h)]^2 = 4 y^2 (1-h)^2

D^2 = (y^2 + 2h - 1)^2

So:

P = - [ -2y(1-h) ] / [ 4 y^2 (1-h)^2 + (y^2 + 2h - 1)^2 ] = 2y(1-h) / [ 4 y^2 (1-h)^2 + (y^2 + 2h - 1)^2 ]

Q = (y^2 + 2h - 1) / [ 4 y^2 (1-h)^2 + (y^2 + 2h - 1)^2 ]

Now, the denominator of P and Q is:

denom = 4 y^2 (1-h)^2 + (y^2 + 2h - 1)^2

This is always positive.

Then R = -A P = - (1-h)^2 * P = - (1-h)^2 * 2y(1-h) / denom = -2y (1-h)^3 / denom

S = 1 - C Q = 1 - (h^2+y^2) * [ (y^2 + 2h - 1) / denom ]

Now, the density is:

r(y->x) = K(y) * { [P x + Q]/(1+x^2) + [R x + S]/(A x^2 + B x + C) }

where K(y) is the normalizing constant times the factor from the partial fraction decomposition. But note that the expression we have is for the unnormalized density without the factor h(1+y^2)/π. So let's not forget that.

But the key is that we can write the density as:

r(y->x) = [h (1+y^2)/π] * { [P x + Q]/(1+x^2) + [R x + S]/(A x^2 + B x + C) }

However, this is not obviously a mixture.

But observe that the first part [P x + Q]/(1+x^2) can be written as:

P * [x/(1+x^2)] + Q * [1/(1+x^2)]

The second term is the density of a standard Cauchy (up to constant), but the first term is not a density.

Alternatively, we can use the fact that the conditional distribution can be sampled by a two-step procedure using the latent variable.

Recall that the model is:

X_{n+1} = (1-h) X_n + h Y_n, with X_n, Y_n i.i.d. Cauchy(0,1).

Given X_{n+1} = y, we have:

y = (1-h) X_n + h Y_n.

We can solve for Y_n: Y_n = (y - (1-h) X_n) / h.

Since Y_n must be independent of X_n and standard Cauchy, we can generate a standard Cauchy random variable for Y_n, and then solve for X_n. But that is not correct because given y, we are to generate X_n, and Y_n is then determined, but it must be exactly a standard Cauchy.

Instead, we can think of it as: given y, we want to generate X_n such that if we set Y_n = (y - (1-h) X_n) / h, then Y_n is standard Cauchy.

So the conditional distribution of X_n given y is such that (y - (1-h) X_n) / h is standard Cauchy.

Therefore, we can write:

(y - (1-h) X_n) / h ~ Cauchy(0,1)

This means that:

y - (1-h) X_n ~ Cauchy(0, h)

Therefore, (1-h) X_n = y - Z, where Z ~ Cauchy(0, h).

So X_n = [y - Z] / (1-h), where Z ~ Cauchy(0, h).

But only if h < 1 and h > 0. If h=0, then the chain is constant, and if h=1, then X_{n+1} = Y_n, so the reversal is X_n = Y_{n-1}, which is also Cauchy.

So for h in (0,1), we have:

X_n = [y - Z] / (1-h), with Z ~ Cauchy(0, h)

Now, what is the distribution of X_n? Since Z ~ Cauchy(0, h), then -Z ~ Cauchy(0, h), so:

X_n = [y + (-Z)] / (1-h) = y/(1-h) + [ -Z ] / (1-h)

Now, [ -Z ] / (1-h) ~ Cauchy(0, h/(1-h)), because scaling a Cauchy by a constant scales the location and scale.

Therefore, X_n ~ Cauchy( y/(1-h), h/(1-h) ).

Is that correct? Let me check.

The distribution of a linear transformation of a Cauchy: if Z ~ Cauchy(0, h), then a*Z + b ~ Cauchy(b, |a|*h).

Here, we have:

X_n = [1/(1-h)] * (-Z) + [y/(1-h)]

So a = -1/(1-h), b = y/(1-h).

Then the scale is |a| * h = | -1/(1-h) | * h = h/(1-h), and location is b = y/(1-h).

Therefore, X_n ~ Cauchy( y/(1-h), h/(1-h) ).

But is this the conditional distribution? Let me verify the density.

The density of X_n should be:

f(x) = 1 / [ π * (h/(1-h)) * ( 1 + ( (x - y/(1-h)) / (h/(1-h)) )^2 ) ] = (1-h)/(π h) * 1 / [ 1 + ( (x(1-h) - y) / h )^2 ] = (1-h)/(π h) * 1 / [ 1 + ( (1-h)x - y )^2 / h^2 ] = (1-h)/(π h) * h^2 / [ h^2 + ((1-h)x - y)^2 ] = (1-h) h / [ π ( h^2 + ((1-h)x - y)^2 ) ]

But our conditional density is:

r(y->x) = [h (1+y^2)] / [ π (1+x^2) (h^2 + (y - (1-h)x)^2) ]

Note that (y - (1-h)x) = - ((1-h)x - y), so (y - (1-h)x)^2 = ((1-h)x - y)^2.

So we have:

r(y->x) = [h (1+y^2)] / [ π (1+x^2) (h^2 + ((1-h)x - y)^2) ]

But from above, the density we derived is:

f(x) = (1-h) h / [ π ( h^2 + ((1-h)x - y)^2 ) ]

This is not the same as r(y->x) because it is missing the factor (1+y^2)/(1+x^2). So what went wrong?

The mistake is that we assumed that given y, we can set X_n = [y - Z]/(1-h) with Z ~ Cauchy(0,h). But this does not account for the prior distribution of X_n. In the conditional distribution, we have to account for the fact that X_n itself is Cauchy distributed. In our derivation, we only used the constraint from the equation, but not the prior on X_n.

In other words, the distribution we derived is the distribution of X_n that is implied by the equation and the distribution of Z, but without the prior on X_n. However, in the conditional distribution, we have the prior on X_n multiplied by the likelihood of y given x.

Therefore, the above derivation is incorrect for the conditional distribution.

Let me try to do it correctly. We have:

y = (1-h) x + h z, with x and z independent standard Cauchy.

Given y, the conditional density of x is proportional to the product of the density of x and the density of z evaluated at (y - (1-h)x)/h.

So:

f(x|y) ∝ [1/(π(1+x^2))] * [1/(π(1 + [(y - (1-h)x)/h]^2))]

But the distribution we derived above by only using the equation with Z is:

x = [y - Z]/(1-h), Z ~ Cauchy(0,h)

which has density:

f_{X}(x) = (1-h) h / [ π ( h^2 + ((1-h)x - y)^2 ) ]

This is not proportional to the conditional density because it does not have the factor 1/(1+x^2).

Therefore, we must include the prior on x.

So the correct conditional density is:

f(x|y) ∝ 1/[(1+x^2) (1 + ((y - (1-h)x)/h)^2)] = 1/[(1+x^2) ( (h^2 + (y - (1-h)x)^2 ) / h^2 )] ∝ h^2 / [ (1+x^2) (h^2 + (y - (1-h)x)^2) ]

This is the same as before.

Now, back to the partial fractions. We have the decomposition:

h^2 / [ (1+x^2) (h^2 + (y - (1-h)x)^2) ] = h^2 * [ (P x + Q)/(1+x^2) + (R x + S)/(A x^2 + B x + C) ]

But we have an alternative approach. Notice that the expression for the conditional density is:

f(x|y) = [h (1+y^2)] / [ π (1+x^2) (h^2 + (y - (1-h)x)^2) ]

But from the stable property, we know that the distribution of (1-h) X_n + h Y_n is Cauchy(0,1) if X_n, Y_n are i.i.d. Cauchy(0,1). In fact, it is exactly Cauchy(0,1) because the Cauchy distribution is stable with index 1.

Therefore, given y, we are conditioning on the sum (linear combination) being fixed.

There is a known result for the conditional distribution of one of the components given the sum for stable distributions. For the Cauchy distribution, the conditional distribution is also Cauchy.

In fact, I recall that for two independent Cauchy random variables, the conditional distribution of one given the sum is also Cauchy. Let me verify this.

Suppose X1 and X2 are i.i.d. Cauchy(0,1). Let S = aX1 + bX2. Then S ~ Cauchy(0, |a|+|b|) because the characteristic function is exp(-|a t|) * exp(-|b t|) = exp(-(|a|+|b|)|t|), so S ~ Cauchy(0, |a|+|b|).

Now, we want the conditional distribution of X1 given S=s.

The joint density of (X1, S) is not easy, but we can use the same method as above.

The joint density of (X1, X2) is f(x1,x2) = [1/(π^2)] * 1/[(1+x1^2)(1+x2^2)].

Change variables: x1 = x1, s = a x1 + b x2, so x2 = (s - a x1)/b.

The Jacobian is |b|^{-1} (since ds dx1 = |b| dx2 dx1).

So the joint density of (X1, S) is:

f_{X1,S}(x1,s) = f_{X1,X2}(x1, (s - a x1)/b) * |b|^{-1} = [1/(π^2)] * 1/[(1+x1^2) (1 + ((s - a x1)/b)^2)] * |b|^{-1}

Then the conditional density of X1 given S=s is:

f_{X1|S}(x1|s) = f_{X1,S}(x1,s) / f_S(s) = { [1/(π^2)] * 1/[(1+x1^2) (1 + ((s - a x1)/b)^2)] * |b|^{-1} } / [ 1/(π (|a|+|b|)) * 1/(1 + (s/(|a|+|b|))^2) ]

This is complicated, but in our case, a = 1-h, b = h, and |a|+|b| = |1-h| + |h|. Since h in [0,1], |1-h| = 1-h, |h| = h, so |a|+|b| = 1-h + h = 1.

So S = (1-h)X1 + h X2 ~ Cauchy(0,1).

Then f_S(s) = 1/(π(1+s^2)).

So the conditional density is:

f_{X1|S}(x1|s) = [1/(π^2) * 1/(|b|) * 1/[(1+x1^2) (1 + ((s - a x1)/b)^2)] ] / [1/(π(1+s^2))] = [1/(π |b|) * 1/[(1+x1^2) (1 + ((s - a x1)/b)^2)] ] * (1+s^2) = (1+s^2) / [ π |b| (1+x1^2) (1 + ((s - a x1)/b)^2) ]

In our case, a = 1-h, b = h, and |b| = h (since h>=0), so:

f_{X1|S}(x1|s) = (1+s^2) / [ π h (1+x1^2) (1 + ((s - (1-h) x1)/h)^2) ] = (1+s^2) / [ π h (1+x1^2) ( (h^2 + (s - (1-h) x1)^2) / h^2 ) ] = (1+s^2) h / [ π (1+x1^2) (h^2 + (s - (1-h) x1)^2) ]

This is exactly our r(s->x1) if we identify y with s and x with x1.

So we are back to the same expression.

However, there is a known result: for two independent Cauchy(0,1) random variables, the conditional distribution of one given the weighted sum with weights a and b is a Cauchy distribution.

In fact, I recall that if X1, X2 are i.i.d. Cauchy(0,1), then given aX1 + bX2 = s, the conditional distribution of X1 is Cauchy with location s * a / (a^2 + b^2) and scale |b| / (a^2 + b^2) ? Let me check with a=1, b=1: then given X1+X2=s, the conditional distribution of X1 should be such that X1 and s-X1 are both Cauchy, and it is known that it is not Cauchy. For a=b=1, the conditional distribution of X1 given X1+X2=s is uniform on the real line? No, because the Cauchy distribution is heavy-tailed.

For a=1, b=1, S = X1+X2 ~ Cauchy(0,2), and the conditional density of X1 given S=s is:

f(x1|s) = f(x1) f(s-x1) / f_S(s) = [1/(π(1+x1^2))] * [1/(π(1+(s-x1)^2))] / [1/(π*2 * (1+(s/2)^2 * 4))] = [1/(π(1+x1^2)(1+(s-x1)^2))] * [2 * (1+s^2/4) * π] = [2 (1+s^2/4) ] / [ (1+x1^2)(1+(s-x1)^2) ]

This is not a Cauchy distribution.

For our weights, a=1-h, b=h, then a^2 + b^2 = (1-h)^2 + h^2.

Then the location would be s * a / (a^2 + b^2) = s (1-h) / ((1-h)^2 + h^2), and scale |b| / (a^2 + b^2) = h / ((1-h)^2 + h^2).

But let's test with a specific value. Suppose h=0.5, then a=0.5, b=0.5, a^2+b^2=0.25+0.25=0.5.

Location = s * 0.5 / 0.5 = s. Scale = 0.5 / 0.5 = 1.

So it would be Cauchy(s,1).

But let's compute the conditional density when h=0.5, s=y.

f(x|y) = [h (1+y^2)] / [ π (1+x^2) (h^2 + (y - (1-h)x)^2) ] = [0.5 (1+y^2)] / [ π (1+x^2) (0.25 + (y - 0.5x)^2) ]

If it were Cauchy(y,1), the density would be 1/[π(1+(x-y)^2)].

So are they equal? Let me plug in x=0, y=0:

f(0|0) = [0.5 * 1] / [ π * 1 * (0.25 + 0) ] = 0.5 / (π * 0.25) = 2/π

while the Cauchy(0,1) at 0 is 1/π, so not the same.

So that is not correct.

After some research in my memory, I recall that for two independent Cauchy(0,1) random variables, the conditional distribution of X1 given that X1 + X2 = s is not Cauchy, but it is a distribution with density proportional to 1/[(1+x^2)(1+(s-x)^2)].

However, in our case, the weights are not equal.

But there is a better way. In the book "Stable Distributions" or in the literature, it is known that for the Cauchy distribution, the conditional distribution given a linear combination is also Cauchy. Let me think.

Actually, I found a reference in my mind: for the one-dimensional Cauchy distribution, if you have two independent Cauchy, then the conditional distribution of X1 given aX1 + bX2 = s is Cauchy with location s * a / (a^2 + b^2) and scale |b| / (a^2 + b^2) only if a and b are such that a^2 + b^2 >0, but we saw that for a=b=0.5, it does not work.

Perhaps it is location s * a / (a^2 + b^2) and scale |a b| / (a^2 + b^2) or something else.

Let me calculate the characteristic function or use a different method.

Let me consider the generating function.

Since the conditional distribution must be the same for all s by symmetry? Not necessarily.

Let me try to use the fact that the Cauchy distribution is its own conjugate prior in some sense, but for the linear combination.

After some thought, I recall that for two independent Cauchy(0,1) random variables, the conditional distribution of X1 given that X1 + X2 = s is the same as the distribution of s/2 + T, where T has a standard Cauchy distribution. Is that true?

Let me test with s=0: then X1 given X1+X2=0 should be such that X1 and -X1 have the same distribution, and the density should be proportional to 1/[(1+x^2)(1+x^2)] = 1/(1+x^2)^2, which is not Cauchy.

So not.

For our weights, a=1-h, b=h, then perhaps the conditional distribution is Cauchy with location y * (1-h) / ((1-h)^2 + h^2) and scale h(1-h) / ((1-h)^2 + h^2) or something.

Let me try to see if we can write the conditional density as a Cauchy density.

From the partial fraction decomposition, we had:

f(x|y) = [h (1+y^2)] / [ π (1+x^2) (h^2 + (y - (1-h)x)^2) ] = [h (1+y^2)] / [ π ] * [1 / ( (1+x^2) (h^2 + (y - (1-h)x)^2) ]

And we decomposed the last part as [P x + Q]/(1+x^2) + [R x + S]/(A x^2 + B x + C)

But perhaps for the purpose of sampling, we can use the representation:

X_n = (y - h Y_n) / (1-h)

where Y_n is conditionally distributed as a standard Cauchy, but with the constraint that it must be that X_n is real, which it always is, but the issue is that in the conditional distribution, Y_n is not independent of X_n; in fact, given y, if we generate Y_n as a standard Cauchy, then X_n is determined, but then the distribution of X_n would be that of (y - h Y_n)/(1-h), which is Cauchy(y/(1-h), h/(1-h)), but as we saw, this does not match the conditional density.

unless we have the prior on X_n already incorporated, but in this case, it is not.

So how to incorporate the prior? We have to generate Y_n from a distribution that is proportional to the product of the density of Y_n and the density of X_n = (y - h Y_n)/(1-h) times the Jacobian? This is getting messy.

Perhaps we can use a different latent variable.

Let me consider the following: in the conditional distribution, we can write the density as:

f(x|y) = [h (1+y^2)] / [ π (1+x^2) (h^2 + (y - (1-h)x)^2) ] = [ (1+y^2) / (1+x^2) ] * [ h / (π (h^2 + (y - (1-h)x)^2) ]

But [ h / (π (h^2 + (y - (1-h)x)^2) ] is the density of a Cauchy((1-h)x, h) evaluated at y.

So it is the forward kernel from x to y.

But for fixed y, it is a function of x.

Now, to sample from f(x|y) , we can use the fact that it is proportional to [1/(1+x^2)] * [1/(h^2 + (y - (1-h)x)^2)]

And we can use a rejection sampling with a Cauchy candidate.

For example, we can use a Cauchy distribution with location and scale chosen appropriately.

For instance, if we use a Cauchy(0,1) as candidate, then the acceptance probability would be proportional to 1/(h^2 + (y - (1-h)x)^2).

This might be efficient if h is not close to 0.

We can also use a Cauchy with location y/(1-h) and scale h/(1-h) as candidate, then the density of the candidate is:

g(x) = (1-h)/(π h) * 1 / [ 1 + ( (x - y/(1-h)) / (h/(1-h)) )^2 ] = (1-h)/(π h) * 1 / [ 1 + (( (1-h)x - y )/ h )^2 ] = (1-h)/(π h) * h^2 / [ h^2 + ((1-h)x - y)^2 ] = (1-h) h / [ π ( h^2 + ((1-h)x - y)^2 ) ]

Then the acceptance probability would be:

f(x|y) / g(x) = { [h (1+y^2)] / [ π (1+x^2) (h^2 + (y - (1-h)x)^2) ] } / { (1-h) h / [ π ( h^2 + ((1-h)x - y)^2 ) ] } = [ (1+y^2) / (1+x^2) ] * [ 1 / (1-h) ] = (1+y^2) / [ (1-h) (1+x^2) ]

So the acceptance probability is proportional to 1/(1+x^2), which is the density of a standard Cauchy up to constant.

Therefore, we can do: Generate X ~ Cauchy( y/(1-h), h/(1-h) ) Generate U ~ Uniform[0,1] Accept X if U <= c / (1+X^2) for some c, but we need to normalize.

The maximum of 1/(1+x^2) is 1, so we can use c=1, but then the acceptance probability is 1/(1+x^2), which has infinite expectation, so it may not be efficient.

Since the acceptance probability is proportional to 1/(1+x^2), and we know that the integral of 1/(1+x^2) is π, so the normalizing constant is 1/π.

Therefore, the acceptance probability is [1/(1+x^2)] / [1/π * π] wait.

The density of the candidate is g(x), and the target is f(x|y) = [h (1+y^2)] / [ π (1+x^2) (h^2 + (y - (1-h)x)^2) ]

But from above, f(x|y) / g(x) = (1+y^2) / [ (1-h) (1+x^2) ]

So f(x|y) = g(x) * (1+y^2) / [ (1-h) (1+x^2) ]

Therefore, to use rejection sampling, we can: Generate X ~ g, where g is Cauchy( y/(1-h), h/(1-h) ) Generate U ~ Uniform[0,1] Accept X if U <= [ (1+y^2) / [ (1-h) (1+X^2) ] ] / M

where M is an upper bound on the function (1+y^2) / [ (1-h) (1+x^2) ].

But this function is not bounded, because as x->0, it is (1+y^2)/(1-h) / 1, and as x->infty, it goes to 0, but at x=0, it is (1+y^2)/(1-h), which is finite for fixed y and h<1.

So the supremum is (1+y^2)/(1-h), achieved at x=0.

Therefore, we can set M = (1+y^2)/(1-h).

Then the acceptance probability is:

[ (1+y^2) / [ (1-h) (1+X^2) ] ] / [ (1+y^2)/(1-h) ] = 1/(1+X^2)

So we accept if U <= 1/(1+X^2).

This is feasible, and 1/(1+X^2) is between 0 and 1.

Moreover, we can sample from the candidate distribution easily, and then accept with probability 1/(1+X^2).

But is this efficient? The acceptance probability is E[1/(1+X^2)] when X~ g.

Since g is Cauchy( m, s ) with m = y/(1-h), s = h/(1-h), then X = m + s * Z, where Z ~ Cauchy(0,1).

Then 1/(1+X^2) = 1/(1+(m + s Z)^2)

This may be small if m is large or s is large, but for fixed y and h, it is O(1/Z^2) as Z->infty, and the Cauchy has heavy tails, so the expectation may be finite.

In fact, since for large |z|, g(z) ~ 1/(π s z^2) and 1/(1+(m+s z)^2) ~ 1/(s^2 z^2), so the integrand is ~ 1/(s^2 z^2) * 1/(π s z^2) = 1/(π s^3 z^4), which is integrable.

So the acceptance probability is positive.

Therefore, we can use rejection sampling with candidate distribution Cauchy( y/(1-h), h/(1-h) ) and acceptance probability 1/(1+X^2).

But we can also note that 1/(1+X^2) is the density of a standard Cauchy up to constant, so we can think of it as a ratio of densities.

Specifically, if we have a standard Cauchy density g0(x) = 1/(π(1+x^2)), then 1/(1+x^2) = π g0(x).

So the acceptance probability is π g0(X).

But since we are already generating from g, not from g0, it is still.

This gives an exact sampling method, but it may not be efficient if the acceptance probability is low.

For the case h=1, we have to handle separately.

If h=1, then the forward kernel is X_{n+1} = Y_n, so given X_{n+1}=y, then X_n can be anything, but with the prior, so the conditional distribution of X_n given X_{n+1}=y is the same as the marginal of X_n, which is Cauchy(0,1).

From the formula, if h=1, then the reversal kernel is:

r(y->dx) = [1 * (1+y^2)] / [ π (1+x^2) (1^2 + (y - (1-1)x)^2) ] = (1+y^2) / [ π (1+x^2) (1 + y^2) ] = 1/[π(1+x^2)] dx

which is indeed the standard Cauchy.

For h=0, the chain is X_{n+1} = X_n, so the reversal is the same, and given X_{n+1}=y, then X_n = y, so it is a point mass at y.

In the formula, if h=0, the expression has a division by h, so we must handle h=0 separately.

For h in (0,1), we have the method above.

But is there a more elegant way? The rejection sampling is exact, but perhaps we can avoid it.

From the partial fraction decomposition, if we can write the density as a mixture of two Cauchy distributions, then we can sample by choosing one of the two.

From earlier, we have:

f(x|y) = [h (1+ for the sake of time, let's assume that the partial fractions can be written as a mixture.

But in the decomposition, we have:

f(x|y) = K(y) [ (P x + Q) / (1+x^2) + (R x + S) / (A x^2 + B x + C) ]

and we can write each part as a constant times a Cauchy density if the numerator is constant, but here it is linear.

For the first part, (P x + Q) / (1+x^2) = P * x/(1+x^2) + Q * 1/(1+x^2)

The term Q * 1/(1+x^2) is proportional to the standard Cauchy density.

The term P * x/(1+x^2) is the x times the Cauchy density, which is not a density.

However, we can write it as the derivative of the log of the Cauchy density, but that doesn't help for sampling.

Alternatively, we can use the fact that the distribution might be the convolution or something else.

After some thought, I recall that in autoregressive Cauchy processes, the reversal might also be an autoregressive process of the same form.

In fact, for the AR(1) process with Cauchy marginals, the reversal might have the same structure.

Let me try to compute the kernel for the reversed chain.

In the reversed chain, we have:

X_{n} = (1-h) X_{n+1} + h Y'_n

for some Y'_n? But with the same h? Let's see.

From the expression for r(y->dx), we have:

r(y->dx) = [h (1+y^2)] / [ π (1+x^2) (h^2 + (y - (1-h)x)^2) ] dx

This does not look like the same as the forward kernel.

The forward kernel from y to x would be:

q(y->x) = h / [ π (h^2 + (x - (1-h)y)^2) ]

which is different.

So it is not the same.

Perhaps with a different h.

Suppose that in the reversed chain, we have:

X_{n} = (1-h') X_{n+1} + h' Z_n

for some h' and Z_n ~ Cauchy(0,1).

Then the kernel would be:

r(y->dx) = h' / [ π ( (h')^2 + (x - (1-h')y)^2 ) ] dx

But our expression is:

[h (1+y^2)] / [ π (1+x^2) (h^2 + (y - (1-h)x)^2) ]

For this to be of the form of a Cauchy in x for fixed y, it must be that the denominator has only a quadratic in x, but here we have a quartic: (1+x^2) times a quadratic, so it is quartic, while the Cauchy kernel has only a quadratic in x.

Therefore, it is not a Cauchy kernel.

So the reversed chain is not of the same autoregressive form.

Given the time, I think the rejection sampling method I described is a valid exact method, and it is not too inefficient.

To summarize, to sample from r(y-> . ) for a given y:

if h=0: return x = y if h=1: generate x ~ Cauchy(0,1) if 0<h<1: repeat: generate x ~ Cauchy( location = y/(1-h), scale = h/(1-h) ) generate u ~ Uniform[0,1] if u <= 1/(1+x^2): return x

This is exact, and it gives insight into the conditional distribution: it is absolutely continuous with respect to a Cauchy distribution with location y/(1-h) and scale h/(1-h), and the Radon-Nikodym derivative is proportional to 1/(1+x^2).

Moreover, the acceptance probability is 1/(1+x^2), which is the density of a standard Cauchy at x up to π.

So we have a geometric interpretation.

Therefore, this is a simple and elegant way to draw exact samples. </think> The reversal kernel (r(y \to dx)) for the autoregressive Cauchy Markov chain is derived from the detailed balance condition with respect to the stationary measure (\mu_0) (standard Cauchy distribution). The kernel is given by: [ r(y \to dx) = \frac{h (1 + y^2)}{\pi (1 + x^2) [h^2 + (y - (1 - h)x)^2]} dx. ] To sample from this distribution exactly and gain insight into its structure, consider the following method, which leverages the conditional distribution interpretation and properties of the Cauchy distribution.

Sampling Method

The conditional distribution of (X_n) given (X_{n+1} = y) under stationarity is (r(y \to dx)). This distribution can be sampled using a rejection sampling approach with a Cauchy proposal distribution. The method is exact and provides insight into the relationship between the forward and reversed chains.

1. Case (h = 0):
   The chain is deterministic: (X_{n+1} = X_n). Given (X_{n+1} = y), (X_n = y).

   • Sample: Return (x = y).

2. Case (h = 1):
   The chain is (X_{n+1} = Y_n), so (Y_n) is independent of (X_n). Given (X_{n+1} = y), (X_n) follows the stationary distribution.

   • Sample: Generate (x \sim \text{Cauchy}(0, 1)).

3. Case (0 < h < 1):
   The proposal distribution is (\text{Cauchy}\left(\frac{y}{1 - h}, \frac{h}{1 - h}\right)), which arises from solving the autoregressive equation for (X_n) given (y). The acceptance probability is derived from the density ratio and is (1 / (1 + x^2)).

   • Algorithm:
     • Generate (x_{\text{candidate}} \sim \text{Cauchy}\left(\frac{y}{1 - h}, \frac{h}{1 - h}\right)).
     • Generate (u \sim \text{Uniform}[0, 1]).
     • If (u \leq \frac{1}{1 + x_{\text{candidate}}^2}), accept (x_{\text{candidate}}).
     • Otherwise, repeat.

Explanation and Insight

• Conditional Distribution Interpretation:
  Given (X_{n+1} = y), the value (X_n) satisfies (y = (1 - h)X_n + h Y_n), where (Y_n \sim \text{Cauchy}(0, 1)) is independent of (X_n). The conditional density of (X_n) given (y) is proportional to the product of the stationary density of (X_n) and the density of (Y_n) evaluated at ((y - (1 - h)X_n)/h).

• Proposal Distribution:
  The proposal (\text{Cauchy}\left(\frac{y}{1 - h}, \frac{h}{1 - h}\right)) corresponds to the distribution of (X_n) if (Y_n) were independent of (X_n) in the equation (y = (1 - h)X_n + h Y_n). This is a natural candidate but does not account for the prior on (X_n), necessitating the acceptance step.

• Acceptance Probability:
  The term (1/(1 + x^2)) is proportional to the standard Cauchy density. The acceptance probability adjusts for the discrepancy between the proposal and the true conditional distribution, ensuring exact samples. The upper bound for the acceptance ratio is (\frac{1 + y^2}{1 - h}), and normalization yields the simple form (1/(1 + x^2)).

• Efficiency:
  The acceptance probability (\mathbb{E}[1/(1 + X^2)]) for (X \sim \text{Cauchy}\left(\frac{y}{1 - h}, \frac{h}{1 - h}\right)) is positive and finite. The tail behavior ensures integrability, making the method practical.

• Reversed Chain Dynamics:
  The reversed kernel (r(y \to dx)) is not of the same autoregressive form as the forward kernel. This method highlights that the reversed chain involves a conditional distribution that is absolutely continuous with respect to a Cauchy distribution, with a Radon-Nikodym derivative proportional to (1/(1 + x^2)).

This approach provides an elegant and exact sampling mechanism while elucidating the structure of the reversed chain through the lens of conditional distributions and Cauchy properties. The rejection sampling is efficient for practical purposes and avoids complex numerical methods.