# Probability of the sum of the largest n samples

15 numbers are randomly chosen from U(0,1), what is the probability that the sum of largest four numbers is greater than 3.5?

    With[{f = OrderDistribution[{UniformDistribution[], 15}, #] &},
Probability[a + b + c + d > 3 + 1 / 2,
{a \[Distributed] f, b \[Distributed] f, c \[Distributed] f,
d\[Distributed] f}]]


Trouble is, it is very, very slow. I shut it down after 30 minutes. When I used NProbability it converged on the wrong answer but it did warn me with numerous error messages. Another CAS could do the above code but it also returned the same wrong answer. The right answer is supposed to be (it at least agrees with a simulation):

$\frac{224077804910008595}{584325558976905216}$

$\approx 0.383481094515780$

How do I do this using Mathematica?

• Presumably the distributions of a,b,c and d are linked together. I don't see any such links in your code. To me it looks like you have independent a,b,c,d (from the correct OrderDistributions) when they should be linked. Presume this is why NProbability gives the wrong answer (Probability would give the same but is slower). – Ymareth Jan 5 '14 at 16:14
• @Ymareth I came to the same conclusion, e.g. a must be greater then b. – ybeltukov Jan 5 '14 at 16:17
• @Ymareth and ybeltukov doesn't the fact that it is distributed in f(15) do that? That is the maximum. – bobbym Jan 5 '14 at 16:19
• A possible simulation is num = 10^6; samp = (Sort /@ RandomVariate[UniformDistribution[], {num, 15}])[[All, -4 ;; -1]]; N@Length@Select[samp, Tr@# > 3.5 &]/num The result agrees with the OP's – Dr. belisarius Jan 5 '14 at 17:01
• @belisarius - Agrees with the stated correct answer but how in the OP's code do I know that the f is always greater than f. In your answer you have that implicitly as you're using the joint distribution but in the OP's code there are just 4 independent distributions - unless I'm misreading? – Ymareth Jan 5 '14 at 18:07

## 2 Answers

That's a tricky one!

OrderDistribution[{dist,n}, {k1, k2, ...}] represents the joint (k1, k2,...} th-order statistics distribution from n observations from the distribution dist.

So:

Probability[a + b + c + d > 7/2, {a, b, c, d} \[Distributed]
OrderDistribution[{UniformDistribution[], 15}, {12, 13, 14, 15}]]

(*
224077804910008595/584325558976905216 --> 0.383481
*)

• Hi belisarius; Thanks, can you see why my approach was bad? – bobbym Jan 5 '14 at 17:55
• @bobbym Consider breaking belisarius' code into two independent regions...Probability[a+b+c+d>7/2,{{a,b}[Distributed]OrderDistribution[{UniformDistribution[],15},{12,13}],{c,d}[Distributed]OrderDistribution[{UniformDistribution[],15},{14,15}]}]. This is not the same. Each of these regions can then be split into 4 which will give the same answer as your original code (I think) - run time seems to be very long. – Ymareth Jan 5 '14 at 18:25
• The OPs code is adding four independent random variables. belisarius is instead adding four joint distributions, which is the problem setup described by "the sum of largest four numbers". – bill s Jan 5 '14 at 19:32

A manual approach

Joint (k1, k2,...} th-order statistics distribution is great, but how we can derive the answer by hands?

Let $a, b, c, d$ be the largest 4 number (in any order with each other). Let $g$ be the next largest number. We know that $g\sim \mathop{\rm Beta}(11,1)$. Then the probability under consideration is

NProbability[a + b + c + d > 7/2 \[Conditioned] a > g && b > g && c > g && d > g,
{a \[Distributed] UniformDistribution[], b \[Distributed] UniformDistribution[],
c \[Distributed] UniformDistribution[], d \[Distributed] UniformDistribution[],
g \[Distributed] BetaDistribution[11, 1]}]


0.383481

This is equivalent to the following integral

Binomial[15, 4] NIntegrate[UnitStep[a + b + c + d - 7/2] 11 g^10,
{g, 0, 1}, {a, g, 1}, {b, g, 1}, {c, g, 1}, {d, g, 1}]


0.383481

Binomial[15, 4] comes from the probability that $a, b, c, d$ are the largest 4 number, 11 g^10 is the PDF of the BetaDistribution[11, 1] and the integration limits come from the conditions $a > g$, etc.

Let us consider the indefinite integral

f[a1_, b1_, c1_, d1_] = Integrate[UnitStep[a + b + c + d - 7/2],
{a, -∞, a1}, {b, -∞, b1}, {c, -∞, c1}, {d, -∞, d1}]

1/384 (-7 + 2 a1 + 2 b1 + 2 c1 + 2 d1)^4 UnitStep[-(7/2) + a1 + b1 + c1 + d1]


The definite integral with limits g, 1 in every dimension is (a simple inclusion-exclusion formula)

int[g_] = Total[(-1)^Total /@ Tuples[{0, 1}, 4] f @@@ Tuples[{g, 1}, 4]]

1/384 - 1/96 (-1 + 2 g)^4 UnitStep[-(1/2) + g] +
1/64 (-3 + 4 g)^4 UnitStep[-(3/2) + 2 g] -
1/96 (-5 + 6 g)^4 UnitStep[-(5/2) + 3 g] +
1/384 (-7 + 8 g)^4 UnitStep[-(7/2) + 4 g]


Finally, the probability is

Binomial[15, 4] Integrate[int[g] 11 g^10, {g, 0, 1}]

224077804910008595/584325558976905216