Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

share|improve this question
3  
N@Mean[Max /@ Tuples[Range@6, 2]] –  belisarius Mar 7 '13 at 11:55
3  
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
1  
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
show 2 more comments

3 Answers 3

up vote 19 down vote accepted

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.

share|improve this answer
    
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
4  
@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
add comment

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 *)
share|improve this answer
add comment

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

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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