Michi’s blog » Blog Archive » Response to Heath Raftery

Response to Heath Raftery

February 9th, 2007

While I’m still on the subject of writing code with the PFP library, I may as well join in on a discussion that got pulled into the Carnival of Mathematics exposition.

Heath Raftery writes about weird probabilities in dice discussions, a problem very much reminiscent of the Monty Hall problem, both in the amount of controversy it generates when people discuss it, and also, it seems on at least half the people in the discussion so far, in how much the terms use end up confusing people.

So, in a comment to the Carnival post, the suggestion came from Alex, that a simulation be used to settle the whole question. And since I just wrote some things with PFP, I might as well write another!

The question, as posed by Heath, is

Two dice are rolled simultaneously. Given that one die shows a “4″, what is the probability that the total on the uppermost faces of the two dice is “7″?

so in order to settle it, I’m going to let my computer roll a whole lote of dice. Whenever at least one die shows a 4, I’m going to tally the roll, and whenever the sum is 7, I’m going to tally it as a succesful roll. The quotient of these tallies will be the probabilities.

I end up with the following code

import Probability
import Control.Monad

– A die roll is a uniform random number in [1..6]
die = uniform [1..6::Int]

– Two dice is just a pair of independently rolled die
twoDice = prod die die

tallyOne :: IO (Int, Int)
tallyOne = do
(d1, d2) <- pick twoDice
let (a,b) | d1 == 4 || d2 == 4 = if d1+d2 == 7 then (1,1) else (1,0)
| otherwise = (0,0)
return (a,b)

sumPair (a,b) (c,d) = (a+c,b+d)

simulate n = (fmap (foldr sumPair (0,0)) . sequence . (take n) . repeat) tallyOne

simulation n = fmap (\(x,y) -> (fromIntegral y)/(fromIntegral x)) $ simulate n

which pretty much lets me enter a number of tests to perform, and returns the resulting quotient of tallied successes by tallied valid rolls.

While playing around with it, I discover that at about 100 000 rolls, the execution outgrows the stack, and so I partition the simulation into pieces of 50 000. This yields me the following

*DiceSim> ( fmap (unlines . map show) $ sequence $ take 10
$ repeat $ simulation 50000 ) >>= putStrLn

which then — after several minutes execution time, gives a simulated dataset, namely

0.176424464192368
0.1825110045332107
0.18224936651289714
0.18589282208990382
0.18272166287755764
0.18133159268929505
0.18697201689545934
0.17798213004630536
0.1818062827225131
0.18456002616944717

Plugging it into my OpenOffice calculating tool, I retrieve a mean of 0.1822451368729 and a standard deviation of 0.0032384171164.

Comparing these, now, to the numerical values of
1/6 = 0.166666666666…
2/11 = 0.18181818181818181818…
makes me conclude that sorry, but, the thing that’s within 0.13σ of our experimental mean simply feels much more reliable than that within almost 5σ

Thus, in the end, my bets goes with the teacher of Heath: 2/11, and for about the same reasons as the Monty Hall problem. As an invitation to Heath: if you disagree with my model, please propose how to code or pseudocode the model you interpret, so that we can simulate and test it!

12 People had this to say...

Dan P Said...

I haven’t tested the code but wouldn’t the strict foldl’ in Data.List solve your stack problem with foldr?

February 10th, 2007 at 0:46

Heath Raftery Said...

Love the simulation program, well done.

As I’ve commented back on Alon’s blog, you’ve simulated the alternative question I suggested in my epilogue. Have another read of my original post, and you’ll see that in the last paragraph I describe precisely the scenario you’ve coded – rolls without a 4 are discarded – and agree that the answer to that question is 2/11. My argument is that was not the situation described in the original question, which unfortunately reduces to an argument of ambiguous interpretation. I do enjoy the mathematical argument however – do continue.

February 10th, 2007 at 9:19

Michi Said...

DanP: foldl or even foldl’ make things better. 100 000 tries works without problem.

Heath: I’m putting an answer up at Alon’s.

February 10th, 2007 at 10:48

TBone Said...

Why not just simulate the entire probability space in this case. There are only 6^2=36 possible outcomes given two dice. Each outcome has an associated 1/(6^2)=1/36 probability. That’s a heck of a lot less work than what you are doing and it will give you an exact answer.

You could generate all outcomes pretty quick with list comprehension, and then just fold the sucker into a counting function. Heck, if you use a strict foldl, it won’t even consume any more space than that required to hold one outcome while doing the calculation.

Of course you could just as easily just do the math:

X = Outcome of first dice roll
Y = Outcome of second dice roll

P(X=4&&Y=3 || X=3&&Y=4 given X=4||Y=4) =
P(X=4&&Y=3 || X=3&&Y=4)/P(X=4||Y=4)

P(X=4||Y=4) = P(X=4) + P(Y=4) – P(X=4&&Y=4) = 2/6-1/36 = 11/36
(these events are not mutually exclusive, so you have to subtract one copy of the combined event as otherwise you have counted it twice because it is included in each of the independent events)

P(X=4&&Y=3 || X=3&&Y=4) = P(X=4&&Y=3) + P(X=3&&Y=4) = 2/36
(these two events are mutually exclusive, so nothing has been counted twice — equally, the term you would subtract is P(X=4&&Y=3 && X=3&&Y=4)=0 [being obviously impossible])

Therefore, the answer is 2/36 * 36/11 = 2/11.

February 10th, 2007 at 15:13

Michi Said...

TBone: I’m pretty certain that’s what the probability library does up to the point that I require, with pick, that it involve random selections.

And the reason I do it this way instead of by actual probability analysis is that the core of the matter was a disagreement in the analysis of the problem, where the original author, Heath, claimed that 2/11 was invalid.

But yeah, the purely mathematical reasoning IS cleaner and better.

February 10th, 2007 at 15:23

TBone Said...

Michl: Yeah, I guess the math argument is not of much use to those who are arguing about the math argument. : )

Interesting post though, I should have a look into that probability library. Sounds like it could be usefull/interesting.

For fun, here’s another (slightly less than optimal) proof in the form ghci commands followed by results.

– The conditional probability space is:

let space = filter (\(x,y) -> x==4 || y==4) [(x,y)|x x+y == 7) space
-- events = [(3,4),(4,3)]

– The probability is therefore (all events in the space being equally likely):

let prob = (length events)/(length space)
– prob = 2/11

February 10th, 2007 at 16:46

TBone Said...

Okay, the blog comment system just totally mulched my code.

Ah well, it was nothing amazing, just a simple list comprehension to generate the die roll space, followed by a filter to reduce that to the conditional space, and another to reduce that to the events space.

Have a good one!

February 10th, 2007 at 16:54

pwyll Said...

The mathematical argument can be made just by counting. There must be at least one 4, so this the sample space: [4 1], [4 2], [4 3], [4 4], [4 5], [4 6], [1 4], [2 4], [3 4] [5 4] [6 4]. There are 11 points in the sample space, only two sum to 7, so the answer is 2/11.

The only trick is to not to count [4 4] more than once. If the reason for this is not obvious, think of the rolls as ordered pairs, the dice as red and green, with 1st dice red, 2nd green.

February 11th, 2007 at 4:59

Michi Said...

pwyll: You did notice which side I’m on in the discussion, right? Heath claims the [4 4] should be double counted. I claim it shouldn’t.

February 11th, 2007 at 12:18

Dan P Said...

Related to this probability question is the Sleeping Beauty Problem. It too depends on asking interesting questions about conditional probability. Unlike Monty Hall and its relatives, however, people still debate about what exactly the correct answer is.

February 12th, 2007 at 19:14

Dan P Said...

I had a look at PFP and don’t like it’s handling of fail. Seems to me that you should be able to write lines in the monad like “guard $ x>10″ so that the final outcome is conditioned on x being greater than 10. The ‘pick’ routine should deal with this by renormalising the probability distribution so as to sum to 1. Maybe I haven’t properly thought through the consequences of this and it’s a stupid suggestion. Does it sound reasonable to you?

February 12th, 2007 at 20:28

Steph Said...

I tried a million iteration of the following program with Matlab and got
>> probparadox
ans =
2.0103 1.0965
>> probparadox
ans =
1.9924 1.0867
The first number is /11 while the second is /6.
Without a doubt, 2/11 wins

n4=0;
n7=0;

for i=1:1000000
x=ceil(6*rand(1));
y=ceil(6*rand(1));
if (x==4 & y~=4)
n4=n4+1;
if x+y==7
n7=n7+1;
end
end
if (y==4 & x~=4)
n4=n4+1;
if x+y==7
n7=n7+1;
end
end
if (x==4 & y==4)
n4=n4+1;
end
end

[(n7/n4)*11 (n7/n4)*6]

February 19th, 2007 at 0:45

Want your say?

post navigation

about: Michi is a recent PhD working in homological algebra and applied algebraic topology. This blog is his outlet for texts with some manner of thought put into them. Over at his LiveJournal intimate details and streams of consciousness might be found.
Not all here is mathematics. All here, though, are my personal thoughts and opinions. Please read the about page (linked above) for more details.
This blog uses statcounter.com for logging and traffic analysis. In order to identify return visitors, this site will issue a cookie on viewing the blog.

Travel plans