Sums of normal random variables

Consider a sample of n independent normal rvs. I would like to identify a systematic way of calculating the probability of having the sum of a subset of them larger than the sum of the rest of rvs. An example case: Population of fish. Mean: 10 kg, stdv: 3 kg. I fish five fish (n=5). What is the probability of having two fish weighing more than the rest of the three fish? The steps which can be followed is to calculate the prob for every combination of fish and then use the inclusion exclusion formula for their union. Is there anything smarter? Note: if four fish were considered the probability of having two of them heavier than the other two should be one. How could this be computed immediately? Thanks for the answers.

• You could certainly do simulation. – Peter Flom Apr 11 at 13:15
• @whuber - You give a great answer assuming that we have a specific two in mind (or randomly choose two). My initial pass at reading thought it was asking about if there were any subsets of 2 such that the sum was greater than the remaining (as evidenced by their claim that if there were 4 fish then the probability would be 1) in which case we would want to look at the distribution of the biggest two vs the distribution of the remaining and would have to dive into the order statistics. Simulation suggests in this situation the probability is roughly .464. – Dason Apr 11 at 14:36
• @Dason Thank you for pointing that out: it is a very plausible interpretation and one I had not conceived of. It also explains why Peter was suggesting simulation, because that's a much trickier problem. I think you're correct about order statistics, because we can reframe the problem as asking "what is the chance that the sum of the $k$ largest of $n$ values exceeds the sum of the $n-k$ smallest ones?" Although we can write down the value as an integral, in general it requires numerical evaluation and rapidly gets onerous as $n$ grows. – whuber Apr 11 at 14:40
• @Manos - If the 1st and 3rd summed were larger than the 2nd, 4th, and 5th... then the 1st and 2nd summed would be larger than the 3rd, 4th, and 5th and would also meet your criteria. So in terms of checking if any subsets meet the criteria we only need to check if the top k sum to something larger than the bottom n-k. – Dason Apr 11 at 16:09
• They could. But s as whuber mentions it's not an easy problem. Simulation would get you a result much easier for any specific situation. – Dason Apr 11 at 16:19

Your example suggests that not only are the $$n$$ variables $$X_1,X_2,\ldots,X_n$$ independent, they also have the same Normal distribution. Let its parameters be $$\mu$$ (the mean) and $$\sigma^2$$ (the variance) and suppose the subset consists of $$k$$ of these variables. We might as well index the variables so that $$X_1,\ldots, X_k$$ are this subset.

The question asks to compute the chance that the sum of the first $$k$$ variables equals or exceeds the sum of the rest:

$$p_{n,k}(\mu,\sigma) = \Pr(X_1+\cdots+X_k \ge X_{k+1}+\cdots+X_n ) = \Pr(Y \le 0)$$

where

$$Y = -(X_1+\cdots+X_k) + (X_{k+1}+\cdots+X_n).$$

$$Y$$ is a linear combination of independent Normal variables and therefore has a Normal distribution--but which one? The laws of expectation and variance immediately tell us

$$E[Y] = -k\mu + (n-k)\mu = (n-2k)\mu$$

and

$$\operatorname{Var}(Y) = k \sigma^2 + (n-k)\sigma^2 = n\sigma^2.$$

Therefore $$Z=\frac{Y - (n-2k)\mu}{\sigma\sqrt{n}}$$ has a standard Normal distribution with distribution function $$\Phi,$$ whence the answer is

$$p_{n,k}(\mu,\sigma) = \Pr(Y \le 0) = \Pr\left(Z \le -\frac{(n-2k)\mu}{\sigma\sqrt{n}}\right) = \Phi\left(-\frac{(n-2k)\mu}{\sigma\sqrt{n}}\right).$$

In the question, $$n=5,k=2,\mu=10,$$ and $$\sigma=3,$$ whence

$$p_{5,2}(10,3) = \Phi\left(-\frac{(5-2(2))10}{3\sqrt{10}}\right)\approx 0.0680186.$$

Generalization

Little needs to change in this analysis even when the $$X_i$$ have different normal distributions or are even correlated: you only need to assume they have an $$n$$-variate Normal distribution to assure their linear combination still has a Normal distribution. The calculations are carried out in the same way and result in a similar formula.

Check

A commenter suggested solving this with simulation. Although that wouldn't be a solution, it's a decent way to check a solution quickly. Thus, in R we might establish the inputs of the simulation in some arbitrary way as

n <- 5
k <- 2
mu <- 10
sigma <- 3
n.sim <- 1e6 # Simulation size
set.seed(17) # For reproducible results

and simulate such data and compare the sums with these two lines:

x <- matrix(rnorm(n*n.sim, mu, sigma), ncol=n)
p.hat <- mean(rowSums(x[, 1:k]) >= rowSums(x[, -(1:k)]))

The post-processing consists of finding the fraction of simulated datasets in which one sum exceeds the other and comparing that to the theoretical solution:

se <- sqrt(p.hat * (1-p.hat) / n.sim)
p <- pnorm(-(n-2*k)*mu / (sigma * sqrt(n)))
signif(c(Simulation=p.hat, Theory=p, Z-score=(p.hat-p)/se), 3)

The output in this case is

Simulation     Theory    Z-score
0.0677     0.0680    -1.1900

The agreement is close and the small absolute z-score allows us to attribute the discrepancy to random fluctuations rather than any error in the theoretical derivation.

• We can also assume without loss of generality that $\sigma=1$; intuitively, we can calculate everything in terms of $\frac {\mu}{\sigma}$ – Acccumulation Apr 11 at 17:33
• @Acccumulation That's correct and it's a good way to proceed. Indeed, this fact follows immediately from observing that one can arbitrarily set the unit of measurement so that $\sigma=1$ without changing the problem. I found it convenient not to have to explain this because it didn't appreciably simplify the analysis. – whuber Apr 11 at 20:15