# Sample a truncated discrete Gaussian distribution?

I need to draw random integers from a truncated discrete Gaussian distribution. That is, the probability of integer $i$ is:

$$p_i \propto \exp\left( -\frac{(i-\mu)^2}{2\sigma^2} \right)$$

for $a \le i \le b$, and $p_i = 0$ otherwise. The missing proportionality constant is obtained from the normalization condition, $\sum_i p_i = 1$. Here $\mu$ and $\sigma$ are real parameters.

Is there a way (package, algorithm, trick, anything) to sample this distribution?

• This is clearly on topic under the help center. I have re-opened. – Glen_b Oct 25 '17 at 0:23
• If the region to sample is not low probability tail region, the simplest approach is to sample from an untruncated normal and throw away the samples that are out of bounds. When the valid region is a high probability region this should work fine, otherwise you will be throwing away a lot of samples and algorithmic time performance may be bad. stats.stackexchange.com/questions/238778/… – Cowboy Trader Oct 25 '17 at 6:23
• For the untruncated case you could use accept/reject ("rejection method"); and it turns out to be surprisingly straightforward - if I have this right you can do it by simply rounding a normal with the same mean and variance as the proposal. You need to scale the pmf of that rounded normal so that it's a proper envelope, but I believe the smallest ratio always occurs at round$(\mu)$, so that's straightforward. It would be possible to reduce the rejection rate further by manipulating the proposal slightly but it's unlikely to be worth it. – Glen_b Oct 25 '17 at 8:06
• Since generating a truncated Gaussian is quite simple (generate a uniform $U$ on $(F(a), F(b))$ and take $F^{-1}(U)$), that should combine with a slight modification of the untruncated scheme above (you need adapt the envelope-scaling-factor $c$ -- in a simple way) to produce a highly efficient method. – Glen_b Oct 25 '17 at 8:10
• @Xi'an Yes, a distribution on the integers. What is confusing on the question? – becko Oct 25 '17 at 8:16

I have trouble with seeing where the difficulty stands in simulating from this integer supported distribution. An R code like

probz=dnorm(a:b,mu,sigma)/sum(dnorm(a:b,mu,sigma))


would work for central values of (a,b), while extreme values would require a bit more caution, e.g.

refz=max(dnorm(a:b,mu,sigma,log=TRUE))
probz=dnorm(a:b,mu,sigma)-refz
probz=exp(probz)/sum(exp(probz))


Then generating by

sample(a:b,n,prob=probz,rep=TRUE)


delivers n iid replicas from this Gaussian-like distribution. Note that if a>mu and b>a are large [in orders of sigma] compared with the mean mu, this distribution will be undistinguishable from a point mass at a.

• Regarding your last sentence, if $\sigma$ is not much smaller than 1, then the probability of $a+1$, $a+2$, and so on for at least a few successive integer values will also be significant – becko Oct 25 '17 at 8:36
• @becko: if $\sigma=1$, $\mu=0$, $a=100<b$, $$p(a)/p(a+1)=10^{44}$$But in any case this is irrelevant for the implementation of the code. – Xi'an Oct 25 '17 at 8:40