# How to model rolling a pair of dice 100 times (Monte Carlo)

Like the title says, I'm having trouble trying to write out a function that will find the probability of finding the sum of a pair of dice 100 times to find the probability of each sum of 2,3,4,5...up to 12.

I've started with

counter[n_][0] = 0; n = 100;
Table[{x[i] = Random[Integer, {1, 6}], y[i] = Random[Integer, {1, 6}]


but I'm not sure how to write the If statement part in which the counter will add 1 for every-time the sum of x[i]+y[i]=2,3,4,5, or 6.

• Tally is the function you're looking for Commented Nov 7, 2015 at 4:26
• Also your Table statement doesn't seem complete. Commented Nov 7, 2015 at 4:28
• You might be interested in this. Commented Nov 8, 2015 at 11:17
• What is the relevance of the 100 times? Are you seeking a sampling distribution of some statistic that depends on the sample size? Commented Nov 8, 2015 at 12:50
• @Wolfies - It's just an example of the Monte Carlo problem by taking many samples and estimating the distribution by the 100 samples. Obviously the bigger the sample size the closer we get to the real distribution. Commented Nov 9, 2015 at 1:42

In Mathematica it is natural to approach such a task with list operations and pattern matching.

dice1 = RandomInteger[{1, 6}, 100];
dice2 = RandomInteger[{1, 6}, 100];
Count[dice1 + dice2, 2 | 3 | 4 | 5 | 6]


You seem to be a very new beginner, since you are using x[i] and y[i] as if these are vectors, when they are in fact not, in Mathematica. Mathematica is very different from other softwares, like MATLAB, so you can't just jump right into it without reading a tutorial first. Nevertheless, here is another solution:

Total@Boole[2 <= Total[#] <= 6 & /@ RandomInteger[{1, 6}, {100, 2}]]

• I'm taking a math modeling class, and my teacher recommended doing it that method since most people in my class have never touched mathematica before. I've taken a basic course on mathematica, but it only included simple stuff. Commented Nov 7, 2015 at 5:15
• Oh also, I forgot to mention that this problem was supposed to be a simulation of the Monte Carlo problem, instead using dice. This was perfect so thanks! Commented Nov 7, 2015 at 5:17
• It may be a matter of taste, but I do find With[ { dist = EmpiricalDistribution @ Map[Total] @ RandomVariate[DiscreteUniformDistribution[{1, 6}], {100, 2}] }, Probability[ 2 <= x <= 6, x \[Distributed] dist] ] more general and easier to understand.
– gwr
Commented Nov 8, 2015 at 16:09
• ... following this I would also consider the use of general statistical functions like RandomVariate and EmpiricalDistribution and Probability the more natural approach in Mathematica given what it is? (cf. my answer)
– gwr
Commented Nov 8, 2015 at 16:13
• @gwr Your answer is very nice (you certainly have my vote), although I don't think the OP is trying to learn the probability functions in Mathematica. The way I interpreted the question he is just trying to learn basic Mathematica skills, and I chose my approach based on that. Commented Nov 8, 2015 at 20:57

There are many ways you can do this, e.g.

ri = RandomInteger[{1, 6}, {100, 2}];
SortBy[Normal@GroupBy[ri, Total, Length@#/100. &], First]


yielding:

{2 -> 0.01, 3 -> 0.05, 4 -> 0.11, 5 -> 0.08, 6 -> 0.13, 7 -> 0.12, 8 -> 0.17, 9 -> 0.14, 10 -> 0.11, 11 -> 0.07, 12 -> 0.01}

You can also exploit DiscreteUniformDistribution and TransformedDistribution, e.g.:

dis = DiscreteUniformDistribution[{1, 6}];
td = TransformedDistribution[
x + y, {x \[Distributed] dis, y \[Distributed] dis}];


You can sample from td, e.g. RandomVariate[td, 100] and play with.

You can get probability mass function: PDF[td, x]:

and finally just for fun (with polynomials):

func[k_Integer] := Module[{x = Unique[], rules},
rules = CoefficientRules[Sum[x^j, {j, 1, 6}]^2, x]; ({k} /. rules)/
36]
pmf = ProbabilityDistribution[func[k], {k, 2, 12, 1}]


So,

Show[Histogram[RandomVariate[pmf, 10000], Automatic, "PDF"],
DiscretePlot[PDF[pmf, x], {x, 2, 12},
PlotStyle -> {Red, Thickness[0.02]}]]


{#, PDF[pmf, #]} & /@ Range[2, 12] //
TableForm[#, TableHeadings -> {None, {"S", "P(X=S)"}}] &


• +1 for CoefficientRules. Since you have an appreciation for such methods have you seen (40344)? It's my take on methods found on forums.wolfram.com/mathgroup/archive/2011/Feb/msg00226.html. Commented Dec 4, 2015 at 1:10
• @Mr.Wizard thank you. I had not seen and look forward to reading it. :) Commented Dec 4, 2015 at 3:12

While the answers so far have covered a lot of ground already I have not seen EmpiricalDistribution. I would like to build upon this observation by providing a couple of general considerations that I have found to be useful when doing statistical experiments using Mathematica. What users of Mathematica may take for granted may surprise newcomers: You can stay very close to the true "programming language" of models which is Mathematics.

So while one may of course work with low level functions like RandomInteger or Boole or Tally or Count one misses the flexibility and generality of the statistical framework provided by Mathematica which is easily transferable to lots of other cases.

Working with Probabilities and Distributions in Mathematica

I have come to quite like the general way that working with probabilities and distributions is done in Mathematica and whenever possible I try to stay within that framework.

Oversimplifying a bit one might see four sources for a distribution and accordingly separate Mathematica-functions:

1. Parametric Distributions corresponding to some idealized stochastic model (e.g. DiscreteUniformDistribution)
2. Nonparametric Distributions with the most prominent example being that we have a sample from an experiment (e.g. EmipricalDistribution)
3. Formula Distributions where we use ProbabilityDistribution to generate a distribution from a known PDF or CDF
4. Derived Distributions where the distribution is the result of some transformation of random numbers whose distributions are given (e.g. TransformedDistribution

A nice thing to note about (3) is, that for ProbabilityDistribution proportionality suffices which is nice for Bayesian statistics as Mathematica will take care of normalization if given the option Method -> "Normalize").

If one comes up with a distribution, everything from then on will be standard. Thus we can:

So let us look at the question at hand and see how this works.

Finding a Derived Parametric Distribution for the Sum of Two Dice

Given a standard experiment we can immediately provide a parametric distribution that discribes the true (unbiased) distribution asymptotically:

SeedRandom["REPEATABLE@151108"]; (* make everything repeatable *)
\$PlotTheme = "Detailed";

distSingleThrowOneDie = DiscreteUniformDistribution[ {1, 6} ];


From this we can directly derive the parametric distribution for the event "Sum of two dice after a single throw":

distSumTwoDice = TransformedDistribution[
x + y,
{
x \[Distributed] distSingleThrowOneDie,
y \[Distributed] distSingleThrowOneDie
}
];

plotTheoreticalPDF = DiscretePlot[
Evaluate @ PDF[ distSumTwoDice, x ],
{ x, 1, 12 },
ExtentSize -> 0.5,
PlotMarkers -> "Point"
]


We can now nicely use this distribution for the sum of two dice to calculate probabilities:

N @ Probability[ 2 <= x <= 5, x \[Distributed] distSumTwoDice ]


0.277778

Doing Monte Carlo Simulations for Throwing Two Dice

We can use any distribution to sample from it using RandomVariate. So let us throw two dice one million times:

totalSample = With[
{
sampleSize = 1000000
},

(* model throwing dice, so that each die might be given its own \
distribution if needed *)

Transpose @ {
RandomVariate[ distSingleThrowOneDie, sampleSize ],
RandomVariate[ distSingleThrowOneDie, sampleSize ]
}
];


We can now use this experiment to see how the empirical distribution of the sum of two dice changes with growing sample size:

Manipulate[
Module[
{
partialSample,
distEmpiricalSumOfTwoDice
},
(* take the first n results and calculate the Totals *)
partialSample = Map[Total] @ totalSample[[ ;; n ]];
distEmpiricalSumOfTwoDice = EmpiricalDistribution @ partialSample;

(* compare the plots *)
Show @ {
plotTheoreticalPDF,
DiscretePlot[
Evaluate @ PDF[ distEmpiricalSumOfTwoDice, x ],
{x, 2, 12 },
PlotStyle -> Red,
PlotRange -> {{ 1, 13}, {0, 0.3} },
PlotLegends -> None
]
}
],
{{ n, 10},{ 10, 100, 1000, 10000, 100000, 1000000 }}
]


• The answer may be a bit overdone but I felt that some kind of general purpose guide to doing (simple) statistical experiments with Mathematica might be useful - and the question nicely fits as an example. Mathematica offers a very nice coherent framework for working with probability and statistics. Observe, for example, how the PlotLegends (automatically) show the probabilities for the events.
– gwr
Commented Nov 8, 2015 at 15:46
• Nice answer and learning source for me, not overdone at all,
– eldo
Commented Nov 8, 2015 at 21:07
• Prolly not something a novice will easily follow at first, but s/he should be aware Mathematica is capable at this. Commented Nov 8, 2015 at 22:18

You can use RandomVariate to sample from a DiscreteUniformDistribution and then add up the pairs, calculate the probability of the sums observed, and then extract the probabilities of interest.

(#/100. & /@
Counts[Plus @@@
RandomVariate[
DiscreteUniformDistribution[{1, 6}], {100, 2}]]
)[#] & /@ Range[2, 6]


Hope this helps.

dist = TransformedDistribution[x + y, {
Distributed[x,
DiscreteUniformDistribution[{1, 6}]],
Distributed[y,
DiscreteUniformDistribution[{1, 6}]]}];

SeedRandom[1]


For small sample sizes, the match to the theoretical values is poor.

data = Total /@ RandomInteger[{1, 6}, {100, 2}];

Show[
Histogram[data, {1.5, 12.5, 1}, "PDF",
ChartStyle -> EdgeForm[Gray]],
DiscretePlot[PDF[dist, t], {t, 2, 12},
PlotStyle -> Directive[Red, Thick]]]


The match is good for very large sample sizes..

data = Total /@ RandomInteger[{1, 6}, {10000, 2}];

Show[
Histogram[data, {1.5, 12.5, 1}, "PDF"],
DiscretePlot[PDF[dist, t], {t, 2, 12}]]


Here's a solution for i dice with j faces using IntegerPartitions and Permutations.

dice = 2;
faces = 6;
range = Range[dice, dice*faces];

res1 =
Flatten[Permutations /@ #, 1] & /@
(IntegerPartitions[#, {dice}, Range @ faces] & /@ range);

len = Length /@ res1;

pro1 = 1/faces^dice*len;

Grid @ Join[
{{"Points", "Probability", "Casts", "n Casts"}},
Transpose[{Range[dice, dice*faces], pro1, res1, len}]]


Now we cast the dice many times:

casts = RandomChoice[Flatten[#, 1] &@ res1, 200000];

pro2 = KeySort @ CountsBy[casts, Total]/200000.


<|2 -> 0.027795, 3 -> 0.056105, 4 -> 0.08404, 5 -> 0.111665, 6 -> 0.138945, 7 -> 0.166145, 8 -> 0.137715, 9 -> 0.110435, 10 -> 0.0836, 11 -> 0.05546, 12 -> 0.028095|>

Does this agree with the formula result?

(pro1 - Values @ pro2) // Total


-1.00614*10^-16

Almost!

• A tiny English tip: "dice" is a plural noun, while "die" is singular. Commented Nov 8, 2015 at 9:55
• Oh, thanks, that saves some additional bytes
– eldo
Commented Nov 8, 2015 at 11:12
• Unfortunately "die" is also quite final - maybe not a good choice for repeating something ad infinitum... :)
– gwr
Commented Nov 8, 2015 at 11:24
(* Model an unbiased  six-sided dice throw *)

dice = DiscreteUniformDistribution[{1, 6}];

(* generate random 2 N throws  and compute probabilities of sum *)

Sumdice[j_, num_] :=
Count[RandomVariate[dice, {num, 2}], _?(Total[#] == j &)]/num

ListPlot[Table[{j, N[Sumdice[j, 10^5]]}, {j, 2, 12, 1}],
Filling -> Axis]


In case you want the cumulative probability of getting two different sums (say j and j+2)

Sumdiceb[j_, num_] :=
Count[RandomVariate[
dice, {num, 2}], _?(Total[#] == j \[Or] Total[#] == j + 2 &)]/num

ListPlot[Table[{j, N[Sumdiceb[j, 10^5]]}, {j, 2, 10, 1}],
Filling -> Axis]