Rejection Sampling when the proposal is a mixture of two distributions

Question

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:

Answer 1

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: