14
$\begingroup$

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[15], b \[Distributed] f[14], c \[Distributed] f[13],
      d\[Distributed] f[12]}]]

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?

Improve this question
$\endgroup$
7
  • $\begingroup$ 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). $\endgroup$
    – Ymareth
    Commented Jan 5, 2014 at 16:14
  • $\begingroup$ @Ymareth I came to the same conclusion, e.g. a must be greater then b. $\endgroup$
    – ybeltukov
    Commented Jan 5, 2014 at 16:17
  • $\begingroup$ @Ymareth and ybeltukov doesn't the fact that it is distributed in f(15) do that? That is the maximum. $\endgroup$
    – bobbym
    Commented Jan 5, 2014 at 16:19
  • $\begingroup$ 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 $\endgroup$
    – Dr. belisarius
    Commented Jan 5, 2014 at 17:01
  • $\begingroup$ @belisarius - Agrees with the stated correct answer but how in the OP's code do I know that the f[15] is always greater than f[14]. 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? $\endgroup$
    – Ymareth
    Commented Jan 5, 2014 at 18:07

2 Answers 2

Reset to default
10
$\begingroup$

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
*)
Improve this answer
$\endgroup$
3
  • $\begingroup$ Hi belisarius; Thanks, can you see why my approach was bad? $\endgroup$
    – bobbym
    Commented Jan 5, 2014 at 17:55
  • 1
    $\begingroup$ @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. $\endgroup$
    – Ymareth
    Commented Jan 5, 2014 at 18:25
  • 1
    $\begingroup$ 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". $\endgroup$
    – bill s
    Commented Jan 5, 2014 at 19:32
5
$\begingroup$

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
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.