Question comes from “Introduction to Probability” by Charles M. Grinstead.

Question 4, Section 7.2

Problem Description:

Let \(X\), \(Y\), and \(Z\) be independent random variables with:

\[ \begin{equation} f_X(x) = f_Y(x) = f_Z(x) = \begin{cases} 1 & \text{if } 0 < x < 1\\ 0 & \text{otherwise } \end{cases} \end{equation} \]

Suppose that \(W = X + Y + Z\). Find \(f_W\) directly, and compare your answer with that given by the formula in Example 7.9. Hint: See Example 7.3.

Solution

From example 7.3 we know that if we say \(U = X + Y\) then:

\[ \begin{equation} f_U(u) = \begin{cases} u & \text{if } 0 \leq u \leq 1\\ 2-u & \text{if } 1 \leq u \leq 2 \\ 0 & \text{otherwise} \end{cases} \end{equation} \]

Using this substitution, and a process similar to the one carried out in example 7.2, we can now say that:

\[ \begin{align} f_W(w) &= \int_{-\infty}^{\infty}f_X(w-z-y)f_Y(y) dy \\ f_W(w) &= \int_{-\infty}^{\infty}f_Z(w-u)f_U(u) du \end{align} \]

If \(w - u > 0 \Rightarrow u < w\) or \(w - u < 1 \Rightarrow u > w-1\), then \(f_Z(w-u)=1\), which means that on those bounds:

\[ \begin{align} f_W(w) &= \int_{w-1}^{w}(1)f_U(u) du \\ f_W(w) &= \int_{w-1}^{w}f_U(u) du \end{align} \]

We can now break up our scenario into three cases. Given that the maximum value of \(W\) is three, we check the scenario in which \(0 \leq w \leq 1\), \(1 \leq w \leq 2\), and \(2 \leq w \leq 3\).

Case 1: \(0 \leq w \leq 1\)

If \(0 \leq w \leq 1\), then we know that \(w - 1\) must be less than 0 (and \(w < 1\)). We can then break apart our previous integral into two pieces, representing the part below 0 and the part above:

\[ \begin{align} f_W(w) &= \int_{w-1}^{w}f_U(u) du \\ f_W(w) &= \int_{w-1}^{0}f_U(u) du + \int_{0}^{w}f_U(u) du \end{align} \]

Based our initial definition of \(f_U(u)\), the left integral must be 0, and since we know \(w \leq 1\), the second integral must be equal to \(u\). Thus, we have:

\[ \begin{align} f_W(w) &= 0 + \int_{0}^{w}u du \\ f_W(w) &= \frac{w^2}{2} \end{align} \]

Case 1: \(1 \leq w \leq 2\)

Similarly, if If \(1 \leq w \leq 2\), then we know that \(w - 1\) must be less than 1. Following the same pattern we have:

\[ \begin{align} f_W(w) &= \int_{w-1}^{w}f_U(u) du \\ f_W(w) &= \int_{w-1}^{1}f_U(u) du + \int_{1}^{w}f_U(u) du \end{align} \]

Once again using our initial definition of \(U\), and the fact that \(w \leq 2\) (or :

\[ \begin{align} f_W(w) &= \int_{w-1}^1udu + \int_{1}^{w}(2-u) du \\ f_W(w) &= \frac{u^2}{2}|_{w-1}^1 + 2u-\frac{u^2}{2}|_1^w \\ f_W(w) &= -w^2 + 3w - \frac{3}{2} \end{align} \]

Case 1: \(2 \leq w \leq 3\)

We follow the same steps for this final case:

\[ \begin{align} f_W(w) &= \int_{w-1}^{w}f_U(u) du \\ f_W(w) &= \int_{w-1}^{2}f_U(u) du + \int_{2}^{w}f_U(u) du \\ f_W(w) &= \int_{w-1}^2 (2- u)du + 0 \\ f_W(w) &= 2u - \frac{u^2}{2}|_{w-1}^2 \\ f_W(w) &= \frac{(w-3)^2}{2} \end{align} \]

Putting it all together

Based on the previous three subsections, we have that the final equation for \(f_W(w)\) is:

\[ \begin{equation} f_W(w) = \begin{cases} \frac{w^2}{2} & \text{if } 0 \leq w < 1\\ -w^2 + 3w - \frac{3}{2} & \text{if } 1 \leq w < 2 \\ \frac{(w-3)^2}{2} & \text{if } 2 \leq w \leq 3 \\ 0 & \text{otherwise} \end{cases} \end{equation} \]

Empirical Comparison

library(ggplot2)

We can simulate this scenario by repeatedly choosing random values between 0 and 1 for three different variables, and summing them together:

set.seed(1234)

# use 30,0001 simulations to match seq values in next cell
sims <- 30001
Ws <- numeric(sims)

for (i in 1:3){
  rvs <- runif(sims)
  Ws <- Ws + rvs
}

The cell below creates data using the theoretical solution from above:

range1x <- seq(0, .9999, 0.0001)
range1y <- (range1x * range1x) / 2
range2x <- seq(1, 1.9999, 0.0001)
range2y <- -(range2x * range2x) + 3 * range2x - (3/2)
range3x <- seq(2, 3, 0.0001)
range3y <- ((range3x - 3) ** 2) / 2

xs <- c(range1x, range2x, range3x)
ys <- c(range1y, range2y, range3y)
plt_data <- data.frame(xs, ys, Ws)

A histogram of the simulated data is plotted against the expected density function in the cell below:

ggplot(data=plt_data) + 
  geom_histogram(
    aes(x=Ws, y=..density..), bins = 100, alpha=0.6, color='green') +
  geom_line(aes(x=xs, y=ys), color='red') +
  xlab('w') +
  ylab('F_W(w)')

As is clear from the plot above, the distribution of the simulated values of \(W\) matches the theoretical density function.