Simulating r.v.'s from a joint density by using rejection sampling in R

Question

I wish to sample variables $v$ and $w$ from the joint density $$(v+w)e^{-\frac{(v+w)^{2}}{2x_{0}}-2\mu v-(\mu -\lambda )w},$$ where $x_0$, $\mu$ and $\lambda$ can be seen as positive constant. Since there is a sum $v+w$, thus I replace it with $r=v+w$, and get $$re^{-\frac{r^{2}}{2x_{0}}-(\mu -\lambda)r-(\mu +\lambda)v},$$ then the problem is reduced to sample $r$ and $v$. I use rejection sampling in R to simulate $r$ firstly, then conditional on $r$, simulate $v$.

The marginal density of $r$ is $$f(r) =re^{-\frac{r^{2}}{2x_{0}}-(\mu -\lambda )r}(1-e^{-(\mu +\lambda)r})\frac{1}{\mu +\lambda }$$

Question: the histogram doesn't fit $f(r)$.

I set $Gamma(scale = 1.8, shape = 1.6)$ as a suitable enveloping distribution, $\mu=2$, $\lambda=1$, $x_0=2.5$, the plot and lines of the two densities show that $\frac{f(r)}{Gamma(r,1.8,1.6)}\leqslant\frac{1}{3}$, so I let $M=\frac{1}{3}$. Here is the code.

# Rejection sampling for v and w

# the marginal density of r
mu = 2;  lambda = 1;  x_0 = 2.5

ff <- function(r) {
   r*exp(-r^2/(2*x_0)-(mu-lambda)*r)*(1-exp(-(mu+lambda)*r))/(mu+lambda) }

#   rejection sampling r
r=seq(0,7,0.01)

rejsampling <- function(M) {
while(1) {
u = runif(1, 0, 1)
r = rgamma(1, scale = 1.8, shape = 1.6)

if(dgamma(r,scale = 1.8, shape = 1.6)*u*M <= ff(r))  {
  return(r)
   }
  }
}

M = 0.333
sampr = replicate(10000, rejsampling(M))

hist(sampr, freq=FALSE, ylim=c(0,0.3), main = "f(r) by rejection sampling",breaks=100)
lines(r, ff(r), col = 2)
text(2.5,0.045,labels=paste("f(r)",sep=""),col=2)
lines(r, dgamma(r,scale = 1.8, shape = 1.6))
text(5,0.14,labels="Gamma(scale = 1.8, shape = 1.6)")

But the hist doesn't fit $f(r)$ (the red curve)! Here is the hist. I've checked several times, but couldn't find the problem. So any comment would be the most grateful. Thank you.

Answer 1

I finally find what the problem is, which is the possibility of $f(r)$ being a probability density function. It depends on what parameters, $\mu$, $\lambda$ and $x_0$, we choose. In the adjusted code, I set $\mu=0.5$, $\lambda=1$ and $x_0=0.945$ to assure $f(r)$ being a feasible density, then set the scale factor $M=2$, which work perfectly. Here is the code.

# the marginal density of r
mu = 0.5
lambda = 1
x_0 = 0.945

ff <- function(r) {
r*exp(-r^2/(2*x_0)-(mu-lambda)*r)*(1-exp(-(mu+lambda)*r))/(mu+lambda) }

###########################   rejection sampling r
r=seq(0,7,0.01)
rejsampling <- function(M) {
while(1) {

# Draw single value from between zero and one
u = runif(1, 0, 1)

# Draw candidate value for r from the enveloping distribution
r = rgamma(1, scale = 0.55, shape = 3.5)

# Accept or reject candidate value; if rejected try again
if(dgamma(r,scale = 0.55, shape = 3.5)*u*M <= ff(r))  {
  return(r)
   }
 }
}

M = 2
sampr = replicate(100000, rejsampling(M))

hist(sampr, freq=FALSE, ylim=c(0,0.8), main = "f(r) by rejection sampling",breaks=100,xlab="r")
lines(r, ff(r), col = 2)
text(1.5,0.65,labels=paste("f(r)",sep=""),col=2)
lines(r, dgamma(r,scale = 0.55, shape = 3.5))
text(4,0.25,labels="Gamma(scale=0.55,shape=3.5)")

which gives us a good sampling fit showing as follows,