# Rejection Sampling when the proposal is a mixture of two distributions

Given a density $$f(x; \alpha, 1) = x^{\alpha-1}\exp(-x)/\Gamma(\alpha)$$ and a proposal density $$q(x) = \alpha_1.q_1(x) + \alpha_2.q_2(x)$$ I want to generate random samples using Rejection Sampling method. Here $$q_1(x)$$ and $$q_2(x)$$ are defined as:

$$q_1(x)\propto x^{\alpha-1}.1_{(0,1]}(x)$$

$$q_2(x) \propto \exp(-x).1_{(1, \infty)}(x)$$

The constant $$A$$ for multiplying the proposal density and the probabilities $$\alpha_1$$ and $$\alpha_2$$ for choosing one of the $$q_i$$'s are given.

How should I generate a sample from the $$q(x)$$ which is a mixture of two other densities?

What I have done so far is as follows:

1. I have derived the CDFs of $$q_1$$ and $$q_2$$ and found their inverses.
2. I have generated a U(0,1) sample and chosen one of the $$q_i$$'s based on the values of $$\alpha_1$$ and $$\alpha_2$$ to generate a random sample $$x$$ using Probability Integral Transformation. For that I am generating another $$U(0,1)$$ sample and using the inverse CDF found in the 1st step to generate the sample.
3. Now that I have generated a sample $$x$$ from one of the $$q_i$$'s, I am trying to generate a random sample from $$q(x)$$. For that I am generating another sample $$u$$ from $$U(0,1)$$ and checking whether this $$u <= f(x)/A.q(x)$$ or not. If the condition is satisfied, I am accepting that $$x$$, otherwise rejecting it and repeating the steps until I get an acceptance.

I don't know whether my steps are correct or not. I have tried plotting the true density of a $$Gamma(0.5, 1)$$ and the density generated from my samples, but they look very much different.

Am I doing something wrong? I am using three $$U(0,1)$$ samples. First one is for choosing which $$q_i$$ to use for generating sample, second one is to generate a sample from that $$q_i$$ using PIT and the third one is for generating a sample using the rejection sampling method.

Edited: Here is my R code that I wrote:

set.seed(123)
alpha <- 0.5
alpha1 <- exp(1) / (alpha + exp(1))
alpha2 <- 1 - alpha1
A <- ((1/alpha) + (1/exp(1))) / gamma(alpha)

sample_from_q1 <- function(){
u2 <- runif(1, 0, 1)
return(u2^0.5)
}
sample_from_q2 <- function(){
u2 <- runif(1, 0, 1)
return(1 - log(u2))
}
density_q1 <- function(x){
return ((x>=0 & x<=1) * (alpha * x ^ (alpha - 1)))

}
density_q2 <- function(x){
return ( x>1*(exp(1-x)) )
}
q <- function(x){
qx <- alpha1 * density_q1(x) + alpha2 * density_q2(x)
return (qx)
}
f <- function(x){
return( x^(alpha-1) * exp(-x) / gamma(alpha) )
}

gamma_generator <- function(N){
gamma_samples = rep(NA, N)
for (i in 1:N){
k <- 0
while(TRUE){
u1 <- runif(1, 0, 1)
if(u1<=alpha2){
x <- sample_from_q2()
}
else{
x <- sample_from_q1()
}
u3 <- runif(1, 0, 1)
if(u3 <= f(x)/(A*q(x))){
gamma_samples[i] = x
break
}
k<-k+1
}
}
return(gamma_samples)
}

samples <- gamma_generator(100000)
plot(density((samples)), type = 'l')


The density plot of the generated samples:  • May I suggest once more that you delete the comments as they are not useful for other readers and that you edit your question to provide R plots associated with the current version of the R code rather than R plots produced with an older version. Again, this makes the entry readable by and useful for all others. Jun 5, 2020 at 8:17

The only uncertain issue with your implementation is whether or not you accounted for the normalising constants. Since $$q_1(x)\propto x^{\alpha-1}\mathbb I_{(0,1)}(x)\qquad q_2(x)\propto e^{-x}\mathbb I_{(1,\infty)}(x)$$ we have $$q_1(x)=\alpha x^{\alpha-1}\mathbb I_{(0,1)}(x)\qquad q_2(x)=e^{1-x}\mathbb I_{(1,\infty)}(x)$$ The mixture proposal is then $$q(x)=\alpha_1 \alpha x^{\alpha-1}\mathbb I_{(0,1)}(x) + (1-\alpha_1) e^{1-x}\mathbb I_{(1,\infty)}(x)$$ which can be generated by

1. picking the component (1 versus 2) by generating $$U_1\sim\mathcal U_{(0,1)}$$ and checking whether or not $$U_1<\alpha_1$$
2. generating $$X$$ from either $$q_1$$ as $$X=U_2^{1/\alpha}$$ or $$q_2$$ as $$X=1-\log(U_2)$$ with $$U_2\sim\mathcal U_{(0,1)}$$

The constant $$A$$ is determined by $$q(x)=\alpha_1 \alpha x^{\alpha-1}\mathbb I_{(0,1)}(x) + (1-\alpha_1) e^{1-x}\mathbb I_{(1,\infty)}(x)\ge A^{-1}x^{\alpha-1} e^{-x}/\Gamma(\alpha)$$ Hence $$\frac{\Gamma(a)f(x)}{q(x)}=\begin{cases} e^{-x}/\alpha_1\alpha &\text{ if }x\le 1\\ x^{\alpha-1}/(1-\alpha_1)e &\text{ if }x\ge 1\\ \end{cases}$$ and (imposing $$\alpha\le 1$$) $$A=\max\{ 1/\alpha_1\alpha,1/(1-\alpha_1)e\}\Gamma(\alpha)^{-1}$$ If $$\alpha_1\propto e$$ and $$\alpha_2\propto\alpha$$ then $$\Gamma(\alpha)A=\frac{\alpha+e}{\alpha e}=\frac{1}{\alpha}+\frac{1}{e}$$ Running the accept-reject step means accepting the simulated $$X$$ if $$U_3\le f(X)\big/ A q(X)\qquad U_3\sim\mathcal U_{(0,1)}$$

Checking the method with an R code like

a=0.9 #alpha
A=max(1/a,1/exp(1)) #upper bound A on ratio
gena<-function(T){ #accepted proposals
u=runif(T) #first uniform
y=1-log(u)
x=u^(1/a)
v=runif(T)<.5 #identical weights
x[v]=y[v]
w=runif(T) #third uniform
x[w<x^a/x/exp(x)/((x>1)*exp(1-x)+a*(x<1)*x^a/x)/A]
}
genb<-function(T){ #sample of size T
z=gena(T)
while(length(z)<T)z=c(z,gena(T))
z[1:T]}


gives a good fit to the Gamma$$(.9,1)$$ target density: 