# Sum of normal truncated random variables

Suppose I have $n$ independent normal random variables

$$X_1 \sim \mathrm{N}(\mu_1, \sigma_1^2)\\X_2 \sim \mathrm{N}(\mu_2, \sigma_2^2)\\\vdots\\X_n \sim \mathrm{N}(\mu_n, \sigma_n^2)$$

and $Y=X_1+X_2+\dotsm+X_n$. How would I characterize the density of $Y$ if the distribution of each $X_i$ is each truncated to within $(\mu_i - 2\sigma_i, \mu_i + 2\sigma_i)$? In other words, I'm sampling from $n$ independent normal distributions, discarding samples not within $2\sigma_i$ of each mean, and summing them.

Right now, I'm doing this with the R code below:

x_mu <- c(12, 18, 7)
x_sd <- c(1.5, 2, 0.8)
a <- x_mu - 2 * x_sd
b <- x_mu + 2 * x_sd

samples <- sapply(1:3, function(i) {
return(rtruncnorm(100000, a[i], b[i], x_mu[i], x_sd[i]))
})

y <- rowSums(samples)


Is there any method for generating the density of $Y$ directly?

• Your question implies you know all the $\sigma_i$. Is that really the case or are you estimating them? There's a huge difference! Out of curiosity, why are you throwing away such data? Depending on your objectives, I suspect there exist (much) better procedures. – whuber Jul 3 '14 at 21:14
• I do know all of the means and SDs for my data, yes. – Devin Jul 3 '14 at 21:41
• I believe that you could characterize it as "a mess". This paper, jstor.org/stable/2236545 , examines the matter, with more scientific rigor. – Alecos Papadopoulos Jul 4 '14 at 3:27
• Outside approximation via CLT, this is relatively tricky. I guess if $n$ is small enough you could try numerical convolution. – Glen_b Jul 5 '14 at 7:04
• @Silverfish Depending on implementation, platform and how fine a grid is tolerable, hundreds should be fine (perhaps more); besides speed, though, with enough terms you have to be much more careful about details of implementation or a number of numerical issues can start to crop up. – Glen_b Jan 21 '15 at 15:29

You could use approximation by saddlepoint methods, for the sum of truncated normals. I will not give the details now, you can look at my answer to General sum of Gamma distributions for hints. What we need is to find the moment-generating function for a truncated normal, which is easy. I will do it here for a standard normal truncated at $\pm 2$, which has density $$f(x) =\begin{cases} \frac1{C} \phi(x), & |x| \le 2 \\ 0, & |x| > 2 \end{cases}$$ where $C=\Phi(2) - \Phi(-2)$ here $\phi(x), \Phi(x)$ are density and cdf for a standard normal, respectively.

The moment generating function can be calculated as $$\DeclareMathOperator{\E}{\mathbb{E}} M(t) = \E e^{tX}=\frac1{C}\int_{-2}^2 e^{tx} \phi(x)\; dx=\frac1{C}e^{\frac12 t^2} [\Phi(2-t)-\Phi(-2-t) ]$$ and then we can use saddlepoint approximations.

I'm curious why, but yes, there is a simple way to generate the pdf of this sum of distributions:

## install.packages("truncnorm")
## install.packages("caTools")
library(truncnorm)

x.mu <- c(12, 18, 7)
x.sd <- c(1.5, 2, 0.8)
x.a <- x.mu - 2*x.sd
x.b <- x.mu + 2*x.sd

dmulti <- function(x, a, b, mu, sd)
rowSums(
sapply(1:length(mu),
function(idx)
dtruncnorm(x, a=a[idx], b=b[idx], mean=mu[idx], sd=sd[idx])))/length(mu)
pmulti <- function(q, a, b, mu, sd)
rowSums(
sapply(1:length(mu),
function(idx)
ptruncnorm(q, a=a[idx], b=b[idx], mean=mu[idx], sd=sd[idx])))/length(mu)

pointrange <- range(c(x.a, x.b))
pointseq <- seq(pointrange[1], pointrange[2], length.out=100)
## Plot the probability density function
plot(pointseq, dmulti(pointseq, x.a, x.b, x.mu, x.sd),
type="l")

## Plot the cumulative distribution function
plot(pointseq, pmulti(pointseq, x.a, x.b, x.mu, x.sd),
type="l")

• If I read this code correctly, you seem to be implementing something like a mixture rather than a summation. The plot this code produces is woefully incorrect. It's not even a valid probability density function! – whuber Nov 3 '14 at 18:57
• @whuber, thanks for the catch. I normalized the pdf and added the cdf. – Bill Denney Nov 8 '14 at 12:38
• Thank you. However, the basic error persists: you are computing a mixture distribution rather than the sum. – whuber Nov 8 '14 at 13:17