0
$\begingroup$

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: Density plot of the generated samples

True density generated using dgamma() function in R

$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – Xi'an
    Jun 5, 2020 at 8:17

1 Answer 1

1
$\begingroup$

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:

enter image description here

$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.