What is the average of rolling two dice and only taking the value of the higher dice roll?

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}]]]


and

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$$

-
N@Mean[Max /@ Tuples[Range@6, 2]] –  belisarius Mar 7 '13 at 11:55
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}]]] –  Cameron Murray Mar 7 '13 at 11:56
@CameronMurray: Thank you, this is helpful, yet I really want the analytic solution since I want to play around with it. –  vonjd Mar 7 '13 at 12:20
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 –  Decency Mar 7 '13 at 17:24
@Decency No, I don't think anybody has heard of Python before :-). Did you know this is a Mathematica site? –  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

 Mean[maxDist]


which is the desired 4.47222.

-
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? –  vonjd Mar 7 '13 at 13:22
@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. –  Marcks Thomas 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 *)

-

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}];
die[6]


$\frac{x[1]}{6}+\frac{x[2]}{6}+\frac{x[3]}{6}+\frac{x[4]}{6}+\frac{x[5]}{6}+\frac{x[6]}{6}$

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


$\frac{161}{36}$

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


$\frac{119}{24}$

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

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


$\frac{245}{54}$

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


$\left\{0.0040003,\frac{137}{16}\right\}$

-