# sampling from $\frac{1}{1+x}$ times Gamma distribution density

I am simulating a process by drawing many random variates $$X$$ from a Gamma distribution with parameters $$\alpha$$, $$\beta$$, $$f_X(x) = \frac{\beta^\alpha \, x^{\alpha-1} \, e^{-\beta x}}{\Gamma(\alpha)} \;.$$ The simulation could be made much more efficient if samples $$X$$ were drawn from a modified distribution with p.d.f. $$f_X(x) \propto \frac{1}{1+x} \, \frac{\beta^\alpha \, x^{\alpha-1} \, e^{-\beta x}}{\Gamma(\alpha)}$$ instead. The parameter $$\beta$$ is generally very small so rejection sampling the $$\frac{1}{1+x}$$ factor results in vanishing efficiency. Is there a clever or efficient method for generating the modified samples?

• You might make more progress by asking us about the original question. What is the distribution you really want to sample?
– whuber
Aug 6, 2020 at 13:26
• That distribution is in fact $f_X(x) \propto \frac{1}{1+x} \, \frac{\beta^\alpha \, x^{\alpha-1} \, e^{-\beta x}}{\Gamma(\alpha)}$. The Monte Carlo simulation effectively computes $\int_0^\infty \! \frac{\phi(x)}{1+x} \, \frac{\beta^\alpha \, x^{\alpha-1} \, e^{-\beta x}}{\Gamma(\alpha)} \, dx$ where $\phi(x)$ is some complicated algorithm.
– user293334
Aug 6, 2020 at 13:31
• Okay. Let's try another tack. How general is $\alpha$? Would there be any bounds on it?
– whuber
Aug 6, 2020 at 13:44
• The result would typically be used for inference, where the parameter space for $\alpha$ could reasonably be restricted to the intervals $(0, 1)$ and $(1, 2)$ for two interesting cases.
– user293334
Aug 6, 2020 at 13:50
• You can try to use Metropolis/Metropolis-Hastings Algorithm. Aug 6, 2020 at 17:43

Since $$\frac{1}{1+x} \, \frac{\beta^\alpha \, x^{\alpha-1} \, e^{-\beta x}}{\Gamma(\alpha)}\le \frac{\beta^\alpha \, x^{\alpha-1} \, e^{-\beta x}}{x\Gamma(\alpha)}=\frac{\beta}{\alpha}\frac{\beta^{\alpha-1} \, x^{\alpha-2} \, e^{-\beta x}}{\Gamma(\alpha-1)}$$ another possibility is to accept/reject with a Gamma $$\mathcal G(\alpha-1,\beta)$$ proposal, assuming $$\alpha>1$$. Since the acceptance is driven by the ratio $$x/(1+x)$$, the efficiency would be much improved, compared with $$\frac{1}{1+x} \, \frac{\beta^\alpha \, x^{\alpha-1} \, e^{-\beta x}}{\Gamma(\alpha)}\le \frac{\beta^\alpha \, x^{\alpha-1} \, e^{-\beta x}}{\Gamma(\alpha)}$$ driven by the ratio $$1/(1+x)$$.

For illustration purposes, here is one evaluation of the algorithm:

  e=function(a,b,T=1e5){
for(i in 1:T)
while(runif(1)>1/(1+1/rgamma(1,a,b)))F=F+1
1+F/T}


demonstrating a high acceptance rate for small values of b:

> e(4,.02)
[1] 1.00672
> e(3,.02)
[1] 1.00996
> e(2,.02)
[1] 1.01849
> e(1,.02)
[1] 1.07424
> e(.1,.02)
[1] 3.49856
> e(.1,.001)
[1] 2.14866
> e(.001,.0001)
[1] 116.4172


Note: A much faster check is obtained by computing directly the probability:

f=function(a,b,T=1e6)1/mean(runif(T)<1/(1+1/rgamma(T,a,b)))

• Of course, brilliant! The efficiency of the $1/(1+x)$ sampling is $e^{\beta } \beta ^{\alpha } \Gamma (1-\alpha ,\beta )$, while that of the $x/(1+x)$ sampling is $(\alpha -1) e^{\beta } \beta ^{\alpha -1} \Gamma (1-\alpha ,\beta )$, so the latter is more efficient when $\alpha > 1+\beta$. The remaining problem is then $\alpha$ near unity, where the efficiency goes to zero. (Although your example continues to work?)
– user293334
Aug 7, 2020 at 8:49
• Thanks. It does deteriorate. See my additional run for e(.001,.0001). Aug 7, 2020 at 9:53
• Ah yes, sorry, I didn't understand that a is $\alpha-1$ and not $\alpha$.
– user293334
Aug 7, 2020 at 10:12
• Apologies for the condensed R code! Aug 7, 2020 at 10:34