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

135 Views Asked by At

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 $$f(r,v)=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$, given $0<=v<=r$.

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 of $w$ shows negative value accidentally!

I set $Gamma(scale = 0.55, shape = 3.5)$ as a suitable enveloping distribution, $\mu=0.5$, $\lambda=1$, $x_0=0.945$, $M=2$. 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)")

# Sampling v given r
v = rep(0,N)
x = rep(0,N)

for(i in 1:N) {
   x[i] = rexp(mu+lambda)
   v[i] = x[i]/(1-exp(-(mu+lambda)*sampr[i]))
}

hist(v, freq=FALSE, xlim=c(0,10), ylim=c(0,1), main = "v by rejection sampling",breaks=1000,xlab="v")


##### sampling w
w = rep(0,N)
w = sampr - v
hist(w,100)

Under this algorithm, we have $0<=w<=r$, But, as the following image showing that $w$ has a number of negative values, which should not happen. Any help will be the most grateful. Thank you. enter image description here