# Do the mean and standard deviation of the sum of independent normal random variables $X$ and $Y$ equal $\mu_X+\mu_Y$ and $sd_X+sd_Y$?

I was following with Chapter 4 of famous 'Statistical Inference 2nd ed' textbook (Casella, Berger) and discovered that if X is a randomly normally distributed variable, and so is Y, and both of them are independent, therefore the distribution of the sum of equal length random vectors X + Y has the parameters mu_XY = mu_X + mu_Y and sd_XY = sd_X + sd_Y

Just for check that statement I performed a simulation in R with a following code:

> rnd_vec_X <- rnorm(10000, mean = 10, sd = 3)
> rnd_vec_Y <- rnorm(10000, mean = 20, sd = 7)
> rnd_vec_XY_sum <- rnd_vec_X + rnd_vec_Y
> mean(rnd_vec_XY_sum)
[1] 30.05607
> sd(rnd_vec_XY_sum)
[1] 7.564628


Indeed, the mean of the sum of there independent random variables is ~30 (20 + 10), as predicted by theory. But why the standard deviation is so far from 10 (7 + 3)? I performed nearly twenty simulations with different seeds and different sample sizes... No luck. I guess something is wrong with my code or with my understanding of the topic. But what?

Under independence (in fact, already under lack of correlation), the variance of the sum is the sum of the variances, $$\sigma^2_{X+Y}=\sigma^2_{X}+\sigma^2_{Y}$$ Incidentally, this result does not rely on normality.
It is, therefore, not true that $$\sqrt{\sigma^2_{X+Y}}=\sqrt{\sigma^2_{X}}+\sqrt{\sigma^2_{Y}},$$ i.e., the sum of the standard deviations is not the standard deviation of the sum - simply because the square root of a sum is not the sum of the square roots.
In your example, $X$ and $Y$ have variances 9 and 49, respectively, so that the population standard deviation is
> sqrt(3^2+7^2)