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

Question

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.

Answer 1

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

Answer 2

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}

rules linking sum to frequency.

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)"}}] &

Answer 3

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").

Advantages of Working with Distributions

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

• Calculate Moments
• Use Distribution Related functions like PDF and CDF
• Do Simulation Experiments
• Calculate Probabilities
• Find Expectation Values

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

Answer 4

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.

Answer 5

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

Answer 6

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!

Answer 7

(* 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]