I am trying to sample from a truncated normal using the standard probability integral transform as described on wikipedia with
$$ x = \Phi^{-1} ( \Phi(\alpha) + U*( \Phi(\beta) - \Phi(\alpha)))*\sigma + \mu $$ with $U \sim Uniform(0,1)$, $\alpha = (a-\mu)/\sigma$ and $\beta = (b-\mu)/\sigma$. In my application, $a = 0.0$ and $b = \infty$.
However, this method fails for instance if $\mu = -0.4$ and $\sigma = 0.05$. I know that sampling positive values with these parameters is highly unlikely. Therefore, I guess I run into errors caused by numerical overflow. I use the GSL library in C to sample the values.
Can you recommend any routine to sample values nevertheless? The PIT yields for instance $\infty$ as sample, which I can obviously not process.
This one I already used:
Truncated Normal by John Burkardt
$\begingroup$
$\endgroup$
5
-
$\begingroup$ Sorry, I totally forgot about that post. My bad! I guess I will try the accept-reject method then. I was silently hoping that there exists some kind of library already. For some reason I cannot compile the one by Chopin. $\endgroup$– mscnvrsyCommented Oct 5, 2016 at 15:25
-
$\begingroup$ "Highly unlikely" is somewhat of an understatement; you are trying to sample from a part of the distribution 8 standard deviations from the mean. Even with double precision, to 15 decimal places that probability rounds off to 0.000000000000001. If you draw a million samples a second from a Normal distribution, it would take you about 30 years before you had a 50% chance of drawing one that far (or more) from the mean. $\endgroup$– jbowmanCommented Oct 5, 2016 at 17:16
-
$\begingroup$ Wow, tbh I didn't consider it that way :D I am using the truncated normal as a posterior distribution, which requires a value larger than 0. Would you say it's reasonable to set it equal to the lower boundary (which is 0) in cases where the sampling method fails? $\endgroup$– mscnvrsyCommented Oct 5, 2016 at 19:01
-
$\begingroup$ @Tim Actually I don't think a generalized extreme value distribution is relevant here. $\endgroup$– Glen_bCommented Oct 6, 2016 at 1:15
-
$\begingroup$ @Glen_b you are right. Unfortunately, I can't edit my comment so I deleted it. What I meant to say is that in case of such extreme tails we would rather use extreme value theory in modeling them rather then looking at the ordinary tails. $\endgroup$– TimCommented Oct 6, 2016 at 8:34
Add a comment
|
1 Answer
$\begingroup$
$\endgroup$
1
If you're sampling the extreme tail, a version of the accept reject method X'ian describes in the paper linked in his answer at that earlier question (basically using an exponential proposal) should work really well.
[The comment containing the link was deleted. I mean this earlier Q]
Out this far, the rejection rate should only be about 2.5%; that's pretty decent; you would be generating about 41 exponentials for every 40 truncated normal values you keep.
-
$\begingroup$ Thank you for the illustration! I kind of solved the problem with the negative mean. In a previous publication, the authors decided to block two variables together which led to complicated expressions for the mean. I "unblocked" these variables and get perfectly behaving means now. But thanks to you all for your helpful comments! $\endgroup$– mscnvrsyCommented Oct 6, 2016 at 8:40