I am referencing to this question on mathunderflow:
What is the average of rolling two dice and only taking the value of the higher dice roll?

Now I tried:

Mean[Max[DiscreteUniformDistribution[{1, 6}], DiscreteUniformDistribution[{1, 6}]]]


Mean[Max[DiscreteUniformDistribution[{{1, 6}, {1, 6}}]]]

but both only give me the wrong answer $\frac{7}{2}$.

What shall I do to get the right result:$$E[X] = \sum_{x=1}^6\frac{2(x-1)+1}{36}x = \frac{1}{36}\sum_{x=1}^6(2x^2 - x) = \frac{161}{36} \approx 4.47$$

  • 3
    $\begingroup$ N@Mean[Max /@ Tuples[Range@6, 2]] $\endgroup$ Mar 7 '13 at 11:55
  • 3
    $\begingroup$ I assume you want some kind of simulation of the problem rather than the analytic solution you already have. With a couple of changes to your code N[Mean[Table[ Max[RandomVariate[DiscreteUniformDistribution[{1, 6}]], RandomVariate[DiscreteUniformDistribution[{1, 6}]]], {1000}]]] $\endgroup$ Mar 7 '13 at 11:56
  • $\begingroup$ @CameronMurray: Thank you, this is helpful, yet I really want the analytic solution since I want to play around with it. $\endgroup$
    – vonjd
    Mar 7 '13 at 12:20
  • 1
    $\begingroup$ I thought this was an interesting question so I wrote a quick script in python to calculate the answer, if anyone is familiar with the language: pastebin.com/egvQD3W8 $\endgroup$
    – Decency
    Mar 7 '13 at 17:24
  • 1
    $\begingroup$ @Decency No, I don't think anybody has heard of Python before :-). Did you know this is a Mathematica site? $\endgroup$
    – whuber
    Mar 7 '13 at 20:45

It is possible to do the calculation analytically in Mathematica as well. Here's one way to define the desired distribution:

 maxDist = 
   TransformedDistribution[Max[x, y], 
           {x \[Distributed] DiscreteUniformDistribution[{1, 6}], 
            y \[Distributed] DiscreteUniformDistribution[{1, 6}]}];

Then the mean is calculated


which is the desired 4.47222.

  • $\begingroup$ Thank you, this is what I wanted. Why do you think would my method not work? Or posed differently: When to use TransformedDistribution and when not? $\endgroup$
    – vonjd
    Mar 7 '13 at 13:22
  • 4
    $\begingroup$ @vonjd: Max[a,b] normally remains unevaluated for distributions, but in your method Max[a,a] simplifies to a, regardless of whether a is a number. Effectively, you're taking the mean of one die roll. In this answer, x and y have the same distribution, but are different objects. $\endgroup$ Mar 7 '13 at 15:28

You could create the TransformedDistribution described by bill s or you could use Mathematica's built in OrderDistribution command, which addresses precisely the sort of problem you have posed.

Mean[OrderDistribution[{DiscreteUniformDistribution[{1, 6}], 2}, 2]]//N
(* 4.47222 *)
  • $\begingroup$ Why is k==2 in this case? $\endgroup$
    – Sterling
    Dec 13 '21 at 15:08

For your amusement, here is a cute, Mathematica-oriented solution: represent the die by its probabilities and "teach" the software to combine results using the maximum. So, to begin, here's a way to create dice and an example of the usual six-sided fair die:

Clear[x]; die[n_] = Sum[x[i]/n, {i, n}];


The arguments to x are the possible values and the coefficients in such a sum are their probabilities. Now for the "lessons":

x /: Times[n_x, m_x] := x[Max[Sequence @@ n, Sequence @@ m]];
x /: Power[n_x, m_Integer] := n;
expectation[y_, f_: Identity] := Expand[y] /. x -> f;

The first is the instruction to retain the larger of two values when combining x's via Times (as is conventional); the second states that when multiple x's are combined (via Power) the maximum does not change; and the third is the definition of expectation (when applied to an object like die[6], anyway). The optional argument f is a function whose expectation is requested; by default, it's the expectation of the die's values themselves.

To make these workings evident--they can seem mysterious in practice--let's first show the combination of two dice. It can be done with a multiplication (but keep reading--there are problems with this approach):

die[6] * die[6] // Expand

$\frac{x[1]}{36}+\frac{x[2]}{12}+\frac{5 x[3]}{36}+\frac{7 x[4]}{36}+\frac{x[5]}{4}+\frac{11 x[6]}{36}$

This indicates, for instance, there is a $7/36$ chance that the largest of two independently thrown fair dice is $4$. Finally, the expectation:

die[6] * die[6] // expectation 


By the way, you can use powers to combine multiple identical dice. Want the expectation of the largest when rolling three dice?

die[6]^3 // expectation


What if we were to cap the largest value at $5$?

expectation[die[6]^3, Min[5, #] &]


What about the expected value from the largest of throwing two d10's and a d12?

die[10]^2 * die[12] // expectation

Don't try this! It will take too long to complete. Apparently Expand (as invoked in expectation) suffers from some terribly bad inefficiency. To work around this, use Outer to implement the addition and multiplication that before I had taken for granted:

combine[a_, b_] := Outer[Times, a, b];             (* Multiplies and expands *)
combine[a_, b_, c__] := combine[a, combine[b, c]]; (* Iterates across multiple arguments *)
multiply[a_, n_Integer] /; n >= 1 := Nest[combine[a, #] &, a, n - 1] (* Replaces `Power` *)

Now things are fine:

combine[multiply[die[10], 2], die[12]] // expectation // AbsoluteTiming



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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