I imagine a random process that proceeds as follows:

Suppose there are $n$ dice. ($n$ is pretty big). They exist in an arena. They are all bouncing around. They start in a random configuration. If two dice collide, then the total number of dots after the collision equals the total number of dots before the collision (law of conservation of dots), but apart from this law, the dots are randomized when they collide. After a suitably long time, how many 1's will there be? How many 2's? How many 3's? etc.

The code that I wrote to simulate this process is:

n = 100; 
Dicenumbers = RandomInteger[{1, 6}, n]; 
j = 0; 
While[j < 10000000, 
  g = RandomInteger[{1, n}]; 
  h = RandomInteger[{1, n}];
  Sumnum = 0; 
  Sumnum = Dicenumbers[[h]] + Dicenumbers[[g]]; 
  satisfied = 0; 
  While[satisfied == 0, 
    Dicenumbers[[g]] = RandomInteger[{1, 6}];
    Dicenumbers[[h]] = Sumnum - Dicenumbers[[g]]; 
    If[Dicenumbers[[h]] > 0, If[Dicenumbers[[h]] < 7, satisfied = 1]];];
k = 1; 
m1 = 0; 
m2 = 0;
m3 = 0; 
m4 = 0; 
m5 = 0; 
m6 = 0; 
While[k <= n, 
  If[Dicenumbers[[k]] == 1, m1++]; 
  If[Dicenumbers[[k]] == 2, m2++]; 
  If[Dicenumbers[[k]] == 3, m3++]; 
  If[Dicenumbers[[k]] == 4, m4++]; 
  If[Dicenumbers[[k]] == 5, m5++]; 
  If[Dicenumbers[[k]] == 6, m6++];
BarChart[{m1, m2, m3, m4, m5, m6}]

(I am a pretty bad coder so this is probably really inelegant),

The above code sets the number of dice, $n$, to 100, and creates a bar chart of the dot number distribution after a million random collisions. However, this bar chart only shows the configuration of the dice at the end of the simulation, and I'm more interested into seeing how the number distribution evolves as the collisions go on. My idea was to implement a slider, that showed the state at 0 collisions, then as you scroll it shows the state after $x$ collisions, where $x < 1,000,000$. That was my plan, but I can't figure out how to do it. It doesn't even have to be a bar graph, I'd just be interested in how the system evolves over time.

  • 1
    $\begingroup$ Please post your code as text, and not as a picture, so other people can copy/paste it and test it. See here for a few pointers. $\endgroup$
    – MarcoB
    Aug 7, 2015 at 23:49
  • $\begingroup$ My first concern is that your code does not represent the process you described in words. A simulation of 100 dice should be tracking 600 dots distributed over 100 dice, but your dice seem to have only a single number on them. $\endgroup$
    – m_goldberg
    Aug 8, 2015 at 1:48
  • $\begingroup$ @m_goldberg But the single number on the dice simply represents the number of dots on that dice right? So in effect I still am tracking the 600 dots, just by looking at how many dots are upon each dice(i think). Also - thanks for cleaning up my post! $\endgroup$ Aug 8, 2015 at 1:50
  • $\begingroup$ No, the number of dots on pair of real dice is 30, 15 on each. $\endgroup$
    – m_goldberg
    Aug 8, 2015 at 1:53
  • 1
    $\begingroup$ I would say what you describe in you comment is a much different process than the one you describe in your question. I conclude we should go with your code and forget the informal description, which is more confusing than useful (IMO). $\endgroup$
    – m_goldberg
    Aug 8, 2015 at 2:07

2 Answers 2

(* First we calc the ways to add up n with two dice *)
alts[x_] := Union[Join @@  Map[{#, Reverse@#} &, IntegerPartitions[x, {2}, Range@6]]]
ways = alts /@ Range[1, 12];

(* an initial dice config *)
m1 = RandomVariate[DiscreteUniformDistribution[{1, 6}], 300];

(* the "collision result" function*)
oneCycle[m1_] := Module[{p1, s1, news},
  (* form the collision pairs *)
  p1 = Partition[m1, 2];
  (* calculate the sums to redistribute *)
  s1 = Tr /@ p1;
  (* for each sum chose a way to redistribute *)
  news = RandomChoice /@ ways[[s1]];
  (* shuffle the results for next iteration *)
  k = RandomSample[Flatten@news]];

Monitor[Nest[oneCycle, m1, 1000], Histogram@k]
  • $\begingroup$ I made the probability calcs, and the final distribution is in fact uniform. $\endgroup$ Aug 8, 2015 at 3:24
  • $\begingroup$ But that doesn't seem true for any initial distribution ... $\endgroup$ Aug 8, 2015 at 4:03
  • $\begingroup$ Would Plus be equivalent to Tr here? $\endgroup$
    – BenP1192
    Aug 9, 2015 at 19:39

This is not an answer, but an extended comment on @beliarius's answer which I much admire.

Notwithstanding my admiration, I have some nits to pick.

  1. The simulation should be parameterized by the number of steps to run, the number of dice, and the delay between histograms.
  2. oneCycle should not flatten news twice.
  3. For a million steps, as proposed by the OP, the simulation will run more than 70 hours with a pause of .25 seconds at each histogram.

Therefore, I offer this minor rewrite:

alts[x_] := 
  Union[Join @@ Map[{#, Reverse@#} &, IntegerPartitions[x, {2}, Range @ 6]]]
ways = alts /@ Range[1, 12];

oneCycle[m1_] :=
 Module[{p1, s1, news},
  p1 = Partition[m1, 2];
  s1 = Tr /@ p1;
  news = Flatten[RandomChoice /@ ways[[s1]]];
  k = RandomSample[news, Length @ news]]

In simulation, steps is the number of steps to run, n the number of dice, and dt the pause (in seconds) between histograms.

simulation[steps_Integer: 1000, n_Integer: 100, dt_: .5] :=
  Block[{k = {}, delay = dt},
        RandomVariate[DiscreteUniformDistribution[{1, 6}], n], 
      If[Length @ k > 0, Histogram @ k, ""]];]

Would it be better to gather up the values of k produced for later analysis?


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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