I was wondering how anyone would come up with this idea.
You observe, correctly, that the $Y_i$ have exponential distributions. They were easy to generate from a standard uniform number generator. The question could be put like this:
Find a simple way to exploit your ability to generate $(Y_1,Y_2)$ to draw values $Z$ following any continuous distribution with positive support.
Such a distribution is one with a density function proportional to $e^{g(z)}$ for a function $g$ defined on the positive numbers.
The key terms "simple" and "proportional to" suggest trying a rejection sampling method. That leads to the algorithm in the question in the following generalized form:
Generate $(Y_1,Y_2)$ and keep $Y_1$ provided $Y_2 \gt f(Y_1)$ for some function $f$ to be determined.
Although it might feel more natural to reject when $Y_2\le f(Y_1),$ as we will see this equivalent formulation leads to a simple calculation.
The result of this sampling procedure evidently produces values of $Y_1$ conditional on the event $Y_2 \gt f(Y_1).$ To find its distribution function, apply the (elementary) definition of conditional probability to the event $Y_1 \le z$ for an arbitrary positive number $z.$ It states
$$\Pr(Y_1\le z \mid Y_2 \gt f(Y_1)) \ \propto\ \Pr(Y_1\le z\text{ and }Y_2 \gt f(Y_1)).$$
We needn't be concerned about the constant of proportionality because we can work it out at the very end, knowing the result has to evaluate to $1$ as $z\to\infty$ (by the axiom of Total Probability).
Because $(Y_1,Y_2)$ is independent, their joint density is exponential. Thus, assuming $f(z) \ge 0$ for all $z\gt 0,$
$$\begin{aligned}
\int_0^z e^{g(y_1)}\,\mathrm{d}y_1 &= \Pr(Y_1\le z\text{ and }Y_2 \gt f(Y_1)) \\
&\propto \int_0^z e^{-y_1}\int_{f(y_1)}^\infty e^{-y_2}\,\mathrm{d}y_2\mathrm{d}y_1\\
&= \int_0^z e^{-y_1 - f(y_1)}\,\mathrm{d}y_1.
\end{aligned}$$
Equality will hold for all $z$ provided the two integrands are equal. Solving for $f$ gives
$$f(z) = -z - g(z) + C$$
where the number $C$ accounts for the neglected proportionality constant $e^C.$
Consider the Half Normal distribution where $g(z) = -z^2/2.$ We find
$$f(z) = -z + z^2/2 + C = (1 - z)^2/2 + C - 1/2.$$
We will need $C \ge 1/2$ to assure $f(z)$ is nonnegative. Larger values of $C$ work, too, but cause more rejections and are thereby less efficient.
Clearly, step (3) in the question converts any positive variable (like a Half Normal variable) into a variable symmetrically distributed around $0.$
Applications
For this method to succeed, we need $f$ to attain a minimum that is not too negative. This implies the target distribution must not be too heavy-tailed. One example is the generalized Gamma distribution with density proportional to $\exp(-z^3/3)$ on the positive numbers.
Here are histograms based on a million draws of $(Y_1,Y_2)$ for the Half Normal and Generalized Gamma problems. The red curves plot the target densities to demonstrate the correctness of this algorithm. The (empirical) acceptance rates show how efficient it is.
This R code produced these plots.
set.seed(17)
n <- 1e6
Y <- matrix(-log(runif(2*n)), ncol = 2) # Step (1): obtain iid exponential variates
#
# The function `f`. The constant can be any non-negative value, with 0 being the
# most efficient.
#
Dists <- list(`Half Normal` = function(z, C = 0) (1 - z)^2/2 + C,
`Generalized Gamma` = function(z, C = 0) -z + z^3/3 + 2/3 + C)
pars <- par(mfrow = c(1, length(Dists)))
for(D in names(Dists)) {
f <- Dists[[D]]
z <- Y[Y[, 2] > f(Y[, 1]), 1] # Step (2) of the rejection sampling
rate <- length(z) / nrow(Y)
hist(z, freq = FALSE, main = D,
sub = bquote(paste("Acceptance rate is ", .(signif(rate, 2)))))
g <- function(x) exp(-x - f(x))
A <- integrate(g, 0, Inf)$value # The constant of integration
curve(g(x) / A, add = TRUE, col = "Red", lwd = 2)
}
par(pars)
