Calculating efficiency of rejection sampling I have a joint distribution from which I want to simulate $f(x,y) = 2(1-x)(1-y)(1-xy)^{-3}, 0<x,y<1$.
I have worked out that the marginal distributions are both $U(0,1)$, so I want to simulate by first simulating $Y \sim U(0,1)$, then use rejection sampling to simulate $X$ from $f_{X|Y}$. I have used calculus to work out that for a given $y$, $f_{X|Y}$ is bounded above by $2(1-y)$ if $0<y<\frac{1}{3}$ and $f_{X|Y}$ is bounded above by $\frac{8/27}{y(1-y)}$ if $\frac{1}{3}<y<1$, so I use $U(0,1)$ as the proposal distribution and $M = 2(1-y) \text{ or } \frac{8/27}{y(1-y)}$, depending on $y$.
This seems to work, but I have a problem with calculating the efficiency of the sampling when $\frac{1}{3} < y < 1$, theoretically, it ought to be $\frac{3}{2} \int_{\frac{1}{3}}^{1} \frac{y(1-y)}{8/27} dy \approx 0.625$. However, when I implement it in R, the efficiency (only when $\frac{1}{3} < y < 1$) is approximately $0.2$. I don't understand why this is happening.

XY_pdf <- function(x,y) {
  
  return (2*(1-x)*(1-y)*((1-x*y)^(-3)))
}

simulate_XY <- function() {

  Y = (2/3)*runif(1) + (1/3)
  
  M = 0.296296/((Y)*(1-Y))
  
  candidate_X = runif(1)
  uniform_random_number = runif(1)
  
  counter = 1 
  
  while ((M)*uniform_random_number > XY_pdf(candidate_X, Y)) {
  candidate_X = runif(1)
  uniform_random_number = runif(1)
  counter = counter + 1
}
   
   return(list(candidate_X,Y,counter))
}
  
random_XY <- replicate(1000,simulate_XY())

X_values = unlist(random_XY[1,])
Y_values = unlist(random_XY[2,])
number_of_attempts = unlist(random_XY[3,])

my_df = data.frame(X_values, Y_values, number_of_attempts)

sum(my_df$number_of_attempts)
length(my_df$Y_values)


This doesn't exactly answer the question, but it may be helpful to point out that if $u=\frac{y}{1-y}$ through a transformation of variables, then $u|x\sim\beta'(u;1,2,1,(1-x)^{-1})$, where $\beta'$ is the generalized beta prime distribution.
$u|x$ can be sampled by compounding two gamma distributions after sampling $x\sim{U(0,1)}$ (using the inverse scale parameterization for the gamma distribution):
$$r|x\sim\Gamma(2,(1-x)^{-1})$$
$$u|x\sim\Gamma(1,r)$$
From which we obtain samples of $y|x$ by:
$$y=\frac{u}{u+1}$$
Using this method to sample 1M pairs of $(x,y)$ in R:
x <- runif(1e6)
u <- rgamma(1e6, 1, rgamma(1e6, 2, 1/(1 - x)))
y <- u/(u + 1)

