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?
$\begingroup$
$\endgroup$
5
-
1$\begingroup$ You might make more progress by asking us about the original question. What is the distribution you really want to sample? $\endgroup$– whuber ♦Aug 6, 2020 at 13:26
-
$\begingroup$ 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. $\endgroup$– user293334Aug 6, 2020 at 13:31
-
1$\begingroup$ Okay. Let's try another tack. How general is $\alpha$? Would there be any bounds on it? $\endgroup$– whuber ♦Aug 6, 2020 at 13:44
-
1$\begingroup$ 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. $\endgroup$– user293334Aug 6, 2020 at 13:50
-
1$\begingroup$ You can try to use Metropolis/Metropolis-Hastings Algorithm. $\endgroup$– Koval BorisAug 6, 2020 at 17:43
Add a comment
|
1 Answer
$\begingroup$
$\endgroup$
4
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)))
-
1$\begingroup$ 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?) $\endgroup$– user293334Aug 7, 2020 at 8:49
-
$\begingroup$ Thanks. It does deteriorate. See my additional run for e(.001,.0001). $\endgroup$– Xi'anAug 7, 2020 at 9:53
-
$\begingroup$ Ah yes, sorry, I didn't understand that
ais $\alpha-1$ and not $\alpha$. $\endgroup$– user293334Aug 7, 2020 at 10:12 -