Timeline for How to sample from $c^a d^{a-1} / \Gamma(a)$?
Current License: CC BY-SA 3.0
11 events
| when toggle format | what | by | license | comment | |
|---|---|---|---|---|---|
| Aug 1, 2011 at 19:36 | history | edited | Macro | CC BY-SA 3.0 |
added 237 characters in body
|
| Aug 1, 2011 at 19:35 | comment | added | Macro | I did notice that the "standard deviation rule" that normals follow (68% within 1, 95% within 2, 99.7% within 3) did apply. So basically for large $cd$ it's a point mass at the mode. From what you say, the threshold where this occurs before the numerical problems, so this still works. Thanks for the insight | |
| Aug 1, 2011 at 19:30 | comment | added | whuber♦ | For such large values of $c d$, the distribution essentially is Normal. Its standard deviation is proportional to $\exp(-c d/6)$. For $c d \gt 30$ or so, double precision floats cannot distinguish between this and zero relative to the mode near $\exp(c d)$, so you basically have a delta distribution at the mode! (That's insanely easy to sample from :-). | |
| Aug 1, 2011 at 19:19 | comment | added | Macro | In my experiments I found problems when $cd > 1000$ or so. In that case, the mode of the distribution occurs for larger $x$, so you're evaluating the density at, for example, $x = 1000$. I don't understand how you calculated $e^{cd}$ in the formula you wrote. In my experiments I had to exponentiate $x \log(cd)$ which causes a problem when $x = 1000, cd = 1000$. | |
| Aug 1, 2011 at 19:11 | comment | added | whuber♦ | I have had no difficulties in my experiments. However, to avoid overflow I have from the beginning normalized by dividing $f$ by $f(\exp(c d))$. In other words, exponentiate $a \log(c d) - \log(\Gamma(a)) - \exp(c d) \log(c d) + \log(\Gamma(\exp(c d)))$. You might get underflow, but that's (usually) immaterial: for integration, etc., you can treat underflows as zeros without making any appreciable error. (You do need extended precision arithmetic once $c d$ exceeds 30 or so: almost all the probability is concentrated in such a tiny interval you have to have that precision!) | |
| Aug 1, 2011 at 19:04 | comment | added | Macro | That is what I do for the computation - it still doesn't avoid overflow. You can't exponentiate a number greater than around 500 on a computer. That quantity gets much larger than that. I mean "pretty good" comparing it with the rejection sampling the OP mentioned. | |
| Aug 1, 2011 at 18:42 | comment | added | whuber♦ | A minute for 1,000 variates isn't very good: you will wait hours for one good Monte-Carlo simulation. You can go four orders of magnitude faster using rejection sampling. The trick is to reject with a close approximation of $f$ rather than with respect to a uniform distribution. Concerning the calculation: compute $a \log(c d) - \log(\Gamma(a))$ (by computing log Gamma directly, of course), then exponentiate. That avoids overflow. | |
| Aug 1, 2011 at 17:25 | comment | added | Macro | There is a lot of computing, but it doesn't actually take very long - certainly much faster than rejection sampling. The simulation I showed above took less than a minute. The problem is that when $cd$ is large, it still breaks. This is basically because it has to calculate the equivalent of $(cd)^{x}$ for large $x$. Any solution proposed will have that problem though - I'm trying to figure out if there's a way to do this on the log scale and transforming back. | |
| Aug 1, 2011 at 17:21 | comment | added | whuber♦ | The method is correct, but awfully painful! How many function evaluations do you suppose are needed for a single random variate? Thousands? Tens of thousands? | |
| Aug 1, 2011 at 16:49 | history | edited | Macro | CC BY-SA 3.0 |
added 209 characters in body
|
| Aug 1, 2011 at 16:35 | history | answered | Macro | CC BY-SA 3.0 |