<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/the-deadly-board-game-puzzle-efficient-simulation-in-r/"></div></div><div><a data-url="https://www.r-bloggers.com/the-deadly-board-game-puzzle-efficient-simulation-in-r/" data-counturl="https://www.r-bloggers.com/the-deadly-board-game-puzzle-efficient-simulation-in-r/" data-count="horizontal" data-text="The ‘deadly board game’ puzzle: efficient simulation in R" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://varianceexplained.org/r/board-game-simulation/"> Variance Explained</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>Last Friday’s <a target="_blank" rel="nofollow" href="http://fivethirtyeight.com/features/can-you-survive-this-deadly-board-game/">“The Riddler” column</a> on FiveThirtyEight presents an interesting probabilistic puzzle:</p>
<blockquote>
<p>While traveling in the Kingdom of Arbitraria, you are accused of a heinous crime. Arbitraria decides who’s guilty or innocent not through a court system, but a board game. It’s played on a simple board: a track with sequential spaces numbered from 0 to 1,000. The zero space is marked “start,” and your token is placed on it. You are handed a fair six-sided die and three coins. You are allowed to place the coins on three different (nonzero) spaces. Once placed, the coins may not be moved.</p>
</blockquote>
<blockquote>
<p>After placing the three coins, you roll the die and move your token forward the appropriate number of spaces. If, after moving the token, it lands on a space with a coin on it, you are freed. If not, you roll again and continue moving forward. If your token passes all three coins without landing on one, you are executed. On which three spaces should you place the coins to maximize your chances of survival?</p>
</blockquote>
<p>(There are also two variations offered, for which I present solutions afterwards).</p>
<p>Last year I took a look at <a target="_blank" rel="nofollow" href="http://varianceexplained.org/r/boarding-pass-simulation/">the “lost boarding pass” puzzle</a> in R, and I found it a useful example of using simulation to answer probabilistic puzzles. Much like the boarding pass puzzle, this has recursive relationships between each round and the previous choices, making it nontrivial to simulate. And much like the other puzzle, performing simulations can be a way to gain insight towards an exact solution.</p>
<h3>Simpler version: one choice, 50 spaces</h3>
<p>It’s often good to start with a simpler version of a puzzle. Here, let’s consider the case where:</p>
<ul>
<li>We have only <strong>1</strong> coin to place, not 3</li>
<li>There are only 50 spaces, not 1000</li>
</ul>
<p>Well, to sample a random roll we would use the <code>sample</code> function with <code>replace = TRUE</code>, and can turn those into positions with the <code>cumsum()</code> (cumulative sum) function:</p>
<figure>
<pre><code data-lang="r"><span>set.seed</span><span>(</span><span>2016-10-18</span><span>)</span><span>
</span><span>cumsum</span><span>(</span><span>sample</span><span>(</span><span>6</span><span>,</span><span> </span><span>20</span><span>,</span><span> </span><span>replace</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>))</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">##  [1]  4  7 10 12 13 16 18 23 29 30 35 40 45 46 52 58 59 63 66 71</code></pre>
</figure>
<p>In this simulation the player rolled a 4, then 3, then 3, getting to positions 4, 7, and 10. It’s easy to perform this for many rolls and many trials with the <code>replicate</code> function. I’ll sample 50 rolls in each trial.<sup><a target="_blank" rel="nofollow" href="http://varianceexplained.org/r/board-game-simulation/#fn:fifty">1</a></sup></p>
<figure>
<pre><code data-lang="r"><span>set.seed</span><span>(</span><span>2016-10-18</span><span>)</span><span>
</span><span>num_rolls</span><span> </span><span>&lt;-</span><span> </span><span>50</span><span>
</span><span>max_position</span><span> </span><span>&lt;-</span><span> </span><span>50</span><span>
</span><span>trials</span><span> </span><span>&lt;-</span><span> </span><span>200000</span><span>
</span><span># create a 50 x 10 matrix of cumulative positions
# Each row is a turn and each column a trial
</span><span>positions</span><span> </span><span>&lt;-</span><span> </span><span>replicate</span><span>(</span><span>trials</span><span>,</span><span> </span><span>cumsum</span><span>(</span><span>sample</span><span>(</span><span>6</span><span>,</span><span> </span><span>num_rolls</span><span>,</span><span> </span><span>replace</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>)))</span><span>
</span><span>positions</span><span>[</span><span>1</span><span>:</span><span>6</span><span>,</span><span> </span><span>1</span><span>:</span><span>6</span><span>]</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    4    4    3    2    3    1
## [2,]    7    5    7    7    6    2
## [3,]   10   11   11    8    8    4
## [4,]   12   16   14   14   13    6
## [5,]   13   22   16   18   19    7
## [6,]   16   28   19   23   20    8</code></pre>
</figure>
<p>We end up with 50 rows (one for each roll) and 10 columns (one for each trial), each containing a random simulation. For example, the first player ended up in position 4, then 7, then 10, and so on. The next went from 4 to 5 to 11.</p>
<p>We want to place our one coin in the most likely position, which means we just have to count the number of times each space (up to space 50) is visited. We can use the (built in) <code>tabulate</code> function to do so<sup><a target="_blank" rel="nofollow" href="http://varianceexplained.org/r/board-game-simulation/#fn:tabulate">2</a></sup>:</p>
<figure>
<pre><code data-lang="r"><span>count_per_position</span><span> </span><span>&lt;-</span><span> </span><span>tabulate</span><span>(</span><span>positions</span><span>,</span><span> </span><span>max_position</span><span>)</span><span>
</span><span>count_per_position</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">##  [1] 33078 39284 45305 52811 61869 72038 50617 53673 55955 57916 58285
## [12] 58386 55987 56850 57528 57205 57107 57036 56925 57170 57446 57326
## [23] 57316 56875 57156 57101 57167 57046 57282 57183 57592 56723 57148
## [34] 57281 56971 57088 57274 57256 57543 56971 57217 57194 57186 57159
## [45] 56807 57286 57137 57242 57189 57066</code></pre>
</figure>
<p>I’m going to stick these results into a data frame for easy analysis later:</p>
<figure>
<pre><code data-lang="r"><span>library</span><span>(</span><span>dplyr</span><span>)</span><span>
</span><span>position_probs</span><span> </span><span>&lt;-</span><span> </span><span>data_frame</span><span>(</span><span>position</span><span> </span><span>=</span><span> </span><span>seq_len</span><span>(</span><span>max_position</span><span>),</span><span>
</span><span>probability</span><span> </span><span>=</span><span> </span><span>count_per_position</span><span> </span><span>/</span><span> </span><span>trials</span><span>)</span><span>
</span><span>position_probs</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## # A tibble: 50 × 2
##    position probability
##       &lt;int>       &lt;dbl>
## 1         1    0.165390
## 2         2    0.196420
## 3         3    0.226525
## 4         4    0.264055
## 5         5    0.309345
## 6         6    0.360190
## 7         7    0.253085
## 8         8    0.268365
## 9         9    0.279775
## 10       10    0.289580
## # ... with 40 more rows</code></pre>
</figure>
<h4>Interpreting the results</h4>
<p>Let’s start by visualizing the probability of landing on each space. (I’m using <code>theme_fivethirtyeight</code> from <a target="_blank" rel="nofollow" href="https://github.com/jrnold/ggthemes">ggthemes</a> in honor of the source of the puzzle).</p>
<figure>
<pre><code data-lang="r"><span>library</span><span>(</span><span>ggplot2</span><span>)</span><span>
</span><span>library</span><span>(</span><span>ggthemes</span><span>)</span><span>
</span><span>theme_set</span><span>(</span><span>theme_fivethirtyeight</span><span>()</span><span> </span><span>+</span><span>
</span><span>theme</span><span>(</span><span>axis.title</span><span> </span><span>=</span><span> </span><span>element_text</span><span>()))</span><span>
</span><span>ggplot</span><span>(</span><span>position_probs</span><span>,</span><span> </span><span>aes</span><span>(</span><span>position</span><span>,</span><span> </span><span>probability</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_line</span><span>()</span></code></pre>
</figure>
<p><img alt="center" src="https://i0.wp.com/varianceexplained.org/figs/2016-10-19-board-game-simulation/position_probs_plot-1.png?w=456" data-recalc-dims="1"></p>
<p>First, we notice that each position gets increasingly more likely up to 6, before diving back down. The pattern repeats itself more mildly up to 12, then basically stabilizes.</p>
<p>That the probabilities climb up to 6 makes sense, because there’s only a 1/6 chance of landing on 1 (you have to get it on your first roll), but there are many ways of ending up on 6 (you could roll it directly, roll 1 then 5, two 3s, three 2s, etc). The drop to 7 then makes sense because there is no longer the 1/6 possibility of landing on it on your first roll. So what would be the best spaces to pick if you had only one coin to place?</p>
<figure>
<pre><code data-lang="r"><span>position_probs</span><span> </span><span>%>%</span><span>
</span><span>arrange</span><span>(</span><span>desc</span><span>(</span><span>probability</span><span>))</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## # A tibble: 50 × 2
##    position probability
##       &lt;int>       &lt;dbl>
## 1         6    0.360190
## 2         5    0.309345
## 3        12    0.291930
## 4        11    0.291425
## 5        10    0.289580
## 6        31    0.287960
## 7        39    0.287715
## 8        15    0.287640
## 9        21    0.287230
## 10       22    0.286630
## # ... with 40 more rows</code></pre>
</figure>
<p>6, 5, 12, 11, and 10 (the ends of those periodic cycles) were the only spaces that were better than the stabilizing probability that the later positions get stuck in.</p>
<p>What is that stabilizing state? Well, consider that the average roll of a six-sided die (the average of its six faces) is . This means that in the long run, we’d expect the die to move about 3.5 spaces each round- which means it would hit one out of every 3.5 spots. So if it stabilized around any one number, it would make sense for it to be 1/3.5 (2/7), which looks right:</p>
<p><img alt="center" src="https://i0.wp.com/varianceexplained.org/figs/2016-10-19-board-game-simulation/position_probs_6_12_hline-1.png?w=456" data-recalc-dims="1"></p>
<p>In fact, on reflection it’s fairly straightforward to calculate the exact probabilities for each space by defining a recursive relationship between each position and the previous ones. We could notice that:</p>
<ul>
<li>For each position, the probability is the average of the probabilities of the previous 6 positions</li>
<li>We define the position 0 to have a probability of 1 (you always start there), and negative positions to a have a probability of 0.</li>
</ul>
<p>We can execute this simulation with a for loop or a functional approach<sup><a target="_blank" rel="nofollow" href="http://varianceexplained.org/r/board-game-simulation/#fn:reduce">3</a></sup>:</p>
<figure>
<pre><code data-lang="r"><span># with a for loop:
</span><span>exact</span><span> </span><span>&lt;-</span><span> </span><span>1</span><span>
</span><span>for</span><span> </span><span>(</span><span>i</span><span> </span><span>in</span><span> </span><span>seq_len</span><span>(</span><span>50</span><span>))</span><span> </span><span>{</span><span>
</span><span>exact</span><span> </span><span>&lt;-</span><span> </span><span>c</span><span>(</span><span>exact</span><span>,</span><span> </span><span>sum</span><span>(</span><span>tail</span><span>(</span><span>exact</span><span>,</span><span> </span><span>6</span><span>))</span><span> </span><span>/</span><span> </span><span>6</span><span>)</span><span>
</span><span>}</span><span>
</span><span>exact</span><span> </span><span>&lt;-</span><span> </span><span>exact</span><span>[</span><span>-1</span><span>]</span><span>
</span><span># alternative version with purrr::reduce
</span><span>exact</span><span> </span><span>&lt;-</span><span> </span><span>purrr</span><span>::</span><span>reduce</span><span>(</span><span>seq_len</span><span>(</span><span>51</span><span>),</span><span> </span><span>~</span><span> </span><span>c</span><span>(</span><span>.</span><span>,</span><span> </span><span>sum</span><span>(</span><span>tail</span><span>(</span><span>.</span><span>,</span><span> </span><span>6</span><span>))</span><span> </span><span>/</span><span> </span><span>6</span><span>))[</span><span>-1</span><span>]</span><span>
</span><span>exact</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">##  [1] 0.1666667 0.1944444 0.2268519 0.2646605 0.3087706 0.3602323 0.2536044
##  [8] 0.2680940 0.2803689 0.2892885 0.2933931 0.2908302 0.2792632 0.2835397
## [15] 0.2861139 0.2870714 0.2867019 0.2855867 0.2847128 0.2856211 0.2859680
## [22] 0.2859437 0.2857557 0.2855980 0.2855999 0.2857477 0.2857688 0.2857356
## [29] 0.2857010 0.2856918 0.2857075 0.2857254 0.2857217 0.2857138 0.2857102
## [36] 0.2857117 0.2857151 0.2857163 0.2857148 0.2857137 0.2857136 0.2857142
## [43] 0.2857146 0.2857145 0.2857142 0.2857141 0.2857142 0.2857143 0.2857143
## [50] 0.2857143</code></pre>
</figure>
<p>We can compare this (in red) to our simulated values to confirm we got it right:</p>
<p><img alt="center" src="https://i2.wp.com/varianceexplained.org/figs/2016-10-19-board-game-simulation/unnamed-chunk-5-1.png?w=456" data-recalc-dims="1"></p>
<h3>Three coins to place</h3>
<p>We’ll now try the original version of the problem, where we place three coins. Why not just pick the three individually most likely spaces (5, 6, and 12) and call it a day? Because the probability of landing on these spaces is <strong>not independent</strong>. We need to maximize the chance we land on <strong>any</strong> space if we want to survive, and perhaps choosing 5, 6, and 12 includes “redundancies,” where we are likely.</p>
<p>One common necessity in simulations is to rearrange your data into. Let’s rearrange our simulated data <code>positions</code> into a <strong>binary matrix</strong> that I’m calling <code>boards</code>. Here instead of having one row for each turn, there will be one row for each <strong>space on the board</strong>. We’ll then place 1s at all the spaces where you landed. To do this, we’ll use a trick for using a two-column matrix as indices that we originally used in the <a target="_blank" rel="nofollow" href="http://varianceexplained.org/r/boarding-pass-simulation/">boarding pass post</a>.</p>
<figure>
<pre><code data-lang="r"><span># create a two-column matrix of position and trial number
</span><span>position_indices</span><span> </span><span>&lt;-</span><span> </span><span>cbind</span><span>(</span><span>c</span><span>(</span><span>positions</span><span>),</span><span> </span><span>rep</span><span>(</span><span>seq_len</span><span>(</span><span>trials</span><span>),</span><span> </span><span>each</span><span> </span><span>=</span><span> </span><span>num_rolls</span><span>))</span><span>
</span><span>position_indices</span><span> </span><span>&lt;-</span><span> </span><span>position_indices</span><span>[</span><span>position_indices</span><span>[,</span><span> </span><span>1</span><span>]</span><span> </span><span>&lt;=</span><span> </span><span>max_position</span><span>,</span><span> </span><span>]</span><span>
</span><span># create an empty integer matrix, set those positions to 1
</span><span>boards</span><span> </span><span>&lt;-</span><span> </span><span>matrix</span><span>(</span><span>0L</span><span>,</span><span> </span><span>nrow</span><span> </span><span>=</span><span> </span><span>num_rolls</span><span>,</span><span> </span><span>ncol</span><span> </span><span>=</span><span> </span><span>trials</span><span>)</span><span>
</span><span>boards</span><span>[</span><span>position_indices</span><span>]</span><span> </span><span>&lt;-</span><span> </span><span>1L</span><span>
</span><span>dim</span><span>(</span><span>boards</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## [1]     50 200000</code></pre>
</figure>
<figure>
<pre><code data-lang="r"><span>boards</span><span>[</span><span>1</span><span>:</span><span>10</span><span>,</span><span> </span><span>1</span><span>:</span><span>6</span><span>]</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">##       [,1] [,2] [,3] [,4] [,5] [,6]
##  [1,]    0    0    0    0    0    1
##  [2,]    0    0    0    1    0    1
##  [3,]    0    0    1    0    1    0
##  [4,]    1    1    0    0    0    1
##  [5,]    0    1    0    0    0    0
##  [6,]    0    0    0    0    1    1
##  [7,]    1    0    1    1    0    1
##  [8,]    0    0    0    1    1    1
##  [9,]    0    0    0    0    0    1
## [10,]    1    0    0    0    0    0</code></pre>
</figure>
<p>For example, we can see that the first simulation landed on 4, 7, and 10, and the second on positions 4 and 5.</p>
<p>Why construct it this way? Because <em>it makes it easy to evaluate a three-coin strategy.</em> For example, how well would we do if we choose spaces 1, 2, and 3? We can simply ask how many trials we landed in at least one of them:</p>
<figure>
<pre><code data-lang="r"><span># get first three rows of boards </span><span>mean</span><span>(</span><span>colSums</span><span>(</span><span>boards</span><span>[</span><span>1</span><span>:</span><span>3</span><span>,</span><span> </span><span>])</span><span> </span><span>></span><span> </span><span>0</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## [1] 0.500225</code></pre>
</figure>
<p>It looks like we’d win about half the time (which makes sense- we’d have to get a 1-3 as the first roll). What about choosing 4, 5, and 6?</p>
<figure>
<pre><code data-lang="r"><span>mean</span><span>(</span><span>colSums</span><span>(</span><span>boards</span><span>[</span><span>4</span><span>:</span><span>6</span><span>,</span><span> </span><span>])</span><span> </span><span>></span><span> </span><span>0</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## [1] 0.79342</code></pre>
</figure>
<p>We’d do a good deal better, with an almost 80% chance of surviving.</p>
<p>Let’s brute force this by trying all possible three-coin strategies between 1 and 20. (Since the per-space probability, it looks very unlikely that the best strategies will need to go past 20 spaces. A strategy of picking 5, 6 and 20 isn’t going to be that different than picking 5, 6, and 500). The <a target="_blank" rel="nofollow" href="https://stat.ethz.ch/R-manual/R-devel/library/utils/html/combn.html">combn</a> function is useful for this, creating a matrix with one row for each of the 1140 three-coin choices:</p>
<figure>
<pre><code data-lang="r"><span>strategies</span><span> </span><span>&lt;-</span><span> </span><span>combn</span><span>(</span><span>20</span><span>,</span><span> </span><span>3</span><span>)</span><span>
</span><span>dim</span><span>(</span><span>strategies</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## [1]    3 1140</code></pre>
</figure>
<figure>
<pre><code data-lang="r"><span>strategies</span><span>[,</span><span> </span><span>1</span><span>:</span><span>6</span><span>]</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1    1    1    1    1    1
## [2,]    2    2    2    2    2    2
## [3,]    3    4    5    6    7    8</code></pre>
</figure>
<p>For each of these strategies (each column in the matrix), we do the same “probability-of-winning” calculation as before.</p>
<figure>
<pre><code data-lang="r"><span>calculate_success</span><span> </span><span>&lt;-</span><span> </span><span>function</span><span>(</span><span>choice</span><span>)</span><span> </span><span>mean</span><span>(</span><span>colSums</span><span>(</span><span>boards</span><span>[</span><span>choice</span><span>,</span><span> </span><span>])</span><span> </span><span>></span><span> </span><span>0</span><span>)</span><span>
</span><span># use apply to perform this function on each column of strategies
</span><span>probabilities</span><span> </span><span>&lt;-</span><span> </span><span>apply</span><span>(</span><span>strategies</span><span>,</span><span> </span><span>2</span><span>,</span><span> </span><span>calculate_success</span><span>)</span></code></pre>
</figure>
<p>This is by far the slowest step of the simulation, and it gets slower still if you consider more than 20 spaces. (See the Appendix for a faster matrix-based method when working with two coins).</p>
<p>Which three-coin strategies maximized your probability of surviving?</p>
<figure>
<pre><code data-lang="r"><span># Restructure the strategies matrix into a three-column data frame
</span><span>strategy_df</span><span> </span><span>&lt;-</span><span> </span><span>as_data_frame</span><span>(</span><span>t</span><span>(</span><span>strategies</span><span>))</span><span>
</span><span># now add in the probabilities and sort:
</span><span>best_options</span><span> </span><span>&lt;-</span><span> </span><span>strategy_df</span><span> </span><span>%>%</span><span>
</span><span>mutate</span><span>(</span><span>probability</span><span> </span><span>=</span><span> </span><span>probabilities</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>arrange</span><span>(</span><span>desc</span><span>(</span><span>probability</span><span>))</span><span>
</span><span>best_options</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## # A tibble: 1,140 × 4
##       V1    V2    V3 probability
##    &lt;int> &lt;int> &lt;int>       &lt;dbl>
## 1      4     5     6    0.793420
## 2      5     6     7    0.759905
## 3      3     5     6    0.755580
## 4      5     6     8    0.755060
## 5      5     6     9    0.744740
## 6      5     6    13    0.735625
## 7      5     6    14    0.731195
## 8      5     6    10    0.730910
## 9      5     6    15    0.728470
## 10     5     6    19    0.727695
## # ... with 1,130 more rows</code></pre>
</figure>
<p>We notice that most of the top strategies include 5 and 6 in them, which were the two best choices for our “choose one” version of the game. However, the version that beat everything else (and not by a small amount) was to choose 4-5-6: and even though 4 on its own was not a particularly high probability space.</p>
<p>We can visualize the 25 best choices to see what else they have in common:</p>
<p><img alt="center" src="https://i0.wp.com/varianceexplained.org/figs/2016-10-19-board-game-simulation/best_options_plot-1.png?w=456" data-recalc-dims="1"></p>
<p>We notice that 5 and 6 dominate most strategies, and when we move away from them we tend to focus on consecutive triplets (like 10-11-12 or 9-10-11).</p>
<p>This makes sense in terms of interdependence among choices. Once you’re already picking 5 and 6, you should pick ones that are less likely to co-occur with those two to maximize your chances of hitting at least one. This includes ones that are immediately adjacent (4 or 7).</p>
<p>To review, here’s all the code we used to solve the puzzle, in one place:</p>
<figure>
<pre><code data-lang="r"><span>set.seed</span><span>(</span><span>2016-10-18</span><span>)</span><span>
</span><span>num_rolls</span><span> </span><span>&lt;-</span><span> </span><span>50</span><span>
</span><span>max_position</span><span> </span><span>&lt;-</span><span> </span><span>50</span><span>
</span><span>trials</span><span> </span><span>&lt;-</span><span> </span><span>200000</span><span>
</span><span>positions</span><span> </span><span>&lt;-</span><span> </span><span>replicate</span><span>(</span><span>trials</span><span>,</span><span> </span><span>cumsum</span><span>(</span><span>sample</span><span>(</span><span>6</span><span>,</span><span> </span><span>num_rolls</span><span>,</span><span> </span><span>replace</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>)))</span><span>
</span><span>position_indices</span><span> </span><span>&lt;-</span><span> </span><span>cbind</span><span>(</span><span>c</span><span>(</span><span>positions</span><span>),</span><span> </span><span>rep</span><span>(</span><span>seq_len</span><span>(</span><span>trials</span><span>),</span><span> </span><span>each</span><span> </span><span>=</span><span> </span><span>num_rolls</span><span>))</span><span>
</span><span>position_indices</span><span> </span><span>&lt;-</span><span> </span><span>position_indices</span><span>[</span><span>position_indices</span><span>[,</span><span> </span><span>1</span><span>]</span><span> </span><span>&lt;=</span><span> </span><span>max_position</span><span>,</span><span> </span><span>]</span><span>
</span><span>boards</span><span> </span><span>&lt;-</span><span> </span><span>matrix</span><span>(</span><span>0L</span><span>,</span><span> </span><span>nrow</span><span> </span><span>=</span><span> </span><span>num_rolls</span><span>,</span><span> </span><span>ncol</span><span> </span><span>=</span><span> </span><span>trials</span><span>)</span><span>
</span><span>boards</span><span>[</span><span>position_indices</span><span>]</span><span> </span><span>&lt;-</span><span> </span><span>1L</span><span>
</span><span>calculate_success</span><span> </span><span>&lt;-</span><span> </span><span>function</span><span>(</span><span>choice</span><span>)</span><span> </span><span>mean</span><span>(</span><span>colSums</span><span>(</span><span>boards</span><span>[</span><span>choice</span><span>,</span><span> </span><span>])</span><span> </span><span>></span><span> </span><span>0</span><span>)</span><span>
</span><span>probabilities</span><span> </span><span>&lt;-</span><span> </span><span>apply</span><span>(</span><span>strategies</span><span>,</span><span> </span><span>2</span><span>,</span><span> </span><span>calculate_success</span><span>)</span></code></pre>
</figure>
<h3>Variations</h3>
<p>The Riddler offered a few variations on the puzzle. One of the advantages of the simulation approach to puzzle-solving is that it can be easy to extract the answers to related questions.</p>
<h4>Can’t pick adjacent spaces</h4>
<blockquote>
<p>Suppose there’s an additional rule that you cannot place the coins on adjacent spaces. What is the ideal placement now?</p>
</blockquote>
<p>This is straightforward with our simulation setup:</p>
<figure>
<pre><code data-lang="r"><span>best_options</span><span> </span><span>%>%</span><span>
</span><span>filter</span><span>(</span><span>V</span><span>2</span><span> </span><span>-</span><span> </span><span>V</span><span>1</span><span> </span><span>></span><span> </span><span>1</span><span>,</span><span>
</span><span>V</span><span>3</span><span> </span><span>-</span><span> </span><span>V</span><span>2</span><span> </span><span>></span><span> </span><span>1</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## # A tibble: 816 × 4
##       V1    V2    V3 probability
##    &lt;int> &lt;int> &lt;int>       &lt;dbl>
## 1      6     8    10    0.714270
## 2      4     6     8    0.711690
## 3      4     6    11    0.704025
## 4      4     6     9    0.701600
## 5      6     9    11    0.701550
## 6      4     6    13    0.700615
## 7      6    10    13    0.699085
## 8      4     6    14    0.697805
## 9      6    13    15    0.697280
## 10     4     6    15    0.696115
## # ... with 806 more rows</code></pre>
</figure>
<p>It looks like the best positions are typically those two or three apart, and that include 6.</p>
<h4>Worst spaces</h4>
<blockquote>
<p>What about the worst squares — where should you place your coins if you’re making a play for martyrdom?</p>
</blockquote>
<p>Simply sort your strategies in ascending rather than descending order:</p>
<figure>
<pre><code data-lang="r"><span>best_options</span><span> </span><span>%>%</span><span>
</span><span>arrange</span><span>(</span><span>probability</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## # A tibble: 1,140 × 4
##       V1    V2    V3 probability
##    &lt;int> &lt;int> &lt;int>       &lt;dbl>
## 1      1     2     7    0.475250
## 2      1     2     8    0.500095
## 3      1     2     3    0.500225
## 4      1     3     7    0.501830
## 5      1     2    13    0.516585
## 6      1     7    13    0.520955
## 7      1     2    14    0.522480
## 8      1     2    19    0.523455
## 9      1     2    18    0.523985
## 10     1     2    17    0.524350
## # ... with 1,130 more rows</code></pre>
</figure>
<p>It looks like both 1-2-3 and 1-2-8 offer about fifty-fifty odds, but that 1-2-7 beats them out for the worst combination.</p>
<h3>What if you need to get all three?</h3>
<p>My own variation- what if you needed to land on all three coins to win?</p>
<figure>
<pre><code data-lang="r"><span># change > 0 to == 3 in our simulation
</span><span>calculate_success_all</span><span> </span><span>&lt;-</span><span> </span><span>function</span><span>(</span><span>choice</span><span>)</span><span> </span><span>mean</span><span>(</span><span>colSums</span><span>(</span><span>boards</span><span>[</span><span>choice</span><span>,</span><span> </span><span>])</span><span> </span><span>==</span><span> </span><span>3</span><span>)</span><span>
</span><span>probability_all</span><span> </span><span>&lt;-</span><span> </span><span>apply</span><span>(</span><span>strategies</span><span>,</span><span> </span><span>2</span><span>,</span><span> </span><span>calculate_success_all</span><span>)</span><span>
</span><span>strategy_df</span><span> </span><span>%>%</span><span>
</span><span>mutate</span><span>(</span><span>probability_all</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>arrange</span><span>(</span><span>desc</span><span>(</span><span>probability_all</span><span>))</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## # A tibble: 1,140 × 4
##       V1    V2    V3 probability_all
##    &lt;int> &lt;int> &lt;int>           &lt;dbl>
## 1      6    12    18        0.047035
## 2      6    12    17        0.040835
## 3      5    11    17        0.039325
## 4      6    11    17        0.038795
## 5      6    12    20        0.035380
## 6      6    14    20        0.035085
## 7      8    14    20        0.034615
## 8      6    12    16        0.034585
## 9      5    11    16        0.034475
## 10     4    10    16        0.034210
## # ... with 1,130 more rows</code></pre>
</figure>
<p>The best way to get all three is to pick 6-12-18. In retrospect this makes sense: if you need to hit every coin, this is equivalent to playing the single-choice version <strong>three times</strong> and needing to win all three. Since 6 is the best choice for each “subgame”, you place your coins 6 apart.</p>
<h3>Lessons</h3>
<p>What can we learn from this example simulation?</p>
<ul>
<li>
<p><strong>Start with a simpler version</strong>: In this case we started with a one-coin version of the puzzle, and explored the results with visualization and mathematical reasoning. This gives you an intuition for the problem and help design the full simulation.</p>
</li>
<li>
<p><strong>Add restrictions to your solution space</strong>. The Riddler offered a board with 1,000 spaces. There are 166 million (<code>choose(1000, 3)</code>) possible three-coin strategies in a board that large, which would have made it almost impossible for us to brute-force. But since we’d already started with a simpler verison and knew that the probabilities stabilized around the 15th space, we figured that we didn’t need to use more than the first 20 spaces to find the best strategy. (Exercise for the reader: what is the best strategy that requires a space beyond the first 20, and how highly is it ranked?)</p>
</li>
<li>
<p><strong>Know your built-in functions in R</strong>. The <code>combn</code> function is an easy way to generate possible coin-placing strategies, and the <code>tabulate</code> function is an extremely efficient way of counting integers in a limited range (far faster than the more commonly used <code>table</code>). The <a target="_blank" rel="nofollow" href="http://adv-r.had.co.nz/Vocabulary.html">Vocabulary chapter of Advanced R</a> gives a pared-down list of built-in functions that’s worth reviewing. (Though amusingly, the list skips <code>tabulate</code>!)</p>
</li>
</ul>
<h3>Appendix: Matrix-based approach for two coins</h3>
<p>If you’re experienced in matrix multiplication and sets, you may spot a fast way we can find pairs of “or” combinations.</p>
<p>The number of events in the union of A and B () can be found with:</p>
<p></p>
<p>(where  means the number of times A and B happen together). Luckily, from our binary matrix it is computationally easy to get both  and  as matrices:</p>
<figure>
<pre><code data-lang="r"><span># 50 x 50 matrix of A + B
</span><span>per_position</span><span> </span><span>&lt;-</span><span> </span><span>rowSums</span><span>(</span><span>boards</span><span>)</span><span>
</span><span>position_plus</span><span> </span><span>&lt;-</span><span> </span><span>outer</span><span>(</span><span>per_position</span><span>,</span><span> </span><span>per_position</span><span>,</span><span> </span><span>&quot;+&quot;</span><span>)</span><span>
</span><span># 50 x 50 matrix of A and B
</span><span>position_and</span><span> </span><span>&lt;-</span><span> </span><span>boards</span><span> </span><span>%*%</span><span> </span><span>t</span><span>(</span><span>boards</span><span>)</span><span>
</span><span># 50 x 50 matrix of counts of A or B
</span><span>position_or</span><span> </span><span>&lt;-</span><span> </span><span>position_plus</span><span> </span><span>-</span><span> </span><span>position_and</span></code></pre>
</figure>
<p>Thanks to optimizations in R’s matrix operations, this is far faster than the <code>apply</code> method we use in the main post. We can examine the best strategies using <code>melt</code> from the reshape2 package:</p>
<figure>
<pre><code data-lang="r"><span>library</span><span>(</span><span>reshape2</span><span>)</span><span>
</span><span>melt</span><span>(</span><span>position_or</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>tbl_df</span><span>()</span><span> </span><span>%>%</span><span>
</span><span>filter</span><span>(</span><span>Var2</span><span> </span><span>></span><span> </span><span>Var1</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>mutate</span><span>(</span><span>probability</span><span> </span><span>=</span><span> </span><span>value</span><span> </span><span>/</span><span> </span><span>trials</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>arrange</span><span>(</span><span>desc</span><span>(</span><span>probability</span><span>))</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">## # A tibble: 1,225 × 4
##     Var1  Var2  value probability
##    &lt;int> &lt;int>  &lt;dbl>       &lt;dbl>
## 1      5     6 123477    0.617385
## 2      4     6 114663    0.573315
## 3      6     8 111676    0.558380
## 4      6     9 111578    0.557890
## 5      6    10 110953    0.554765
## 6      6     7 110841    0.554205
## 7      6    13 109681    0.548405
## 8      6    14 109502    0.547510
## 9      6    15 109249    0.546245
## 10     6    31 108944    0.544720
## # ... with 1,215 more rows</code></pre>
</figure>
<p>We can then see that the best two-coin strategies include 5-6 (which are the two best single-coin strategies), as well as 4-6 and 6-8.</p>
<p>I don’t yet have an effective method for extending this matrix approach to three coins. Can you find one?</p>
<div>
<ol>
<li>
<p>We don’t really need to sample 50 rolls for the 50-space version since the chance of needing all fifty (rolling all 1s) is effectively impossible- I include it only for simplicity. <a target="_blank" rel="nofollow" href="http://varianceexplained.org/r/board-game-simulation/#fnref:fifty">↩</a></p>
</li>
<li>
<p>If you’re familiar with R you may have expected me to do <code>table(positions)</code>: but <code>table</code> is <strong>very slow</strong> for large integer vectors relative to other counting methods, taking about two seconds on this data. (Among other reasons, it converts it to a character vector before counting). <code>count</code> from the dplyr package gets about a 10x improvement, and data.table a greater improvement still. But in the very special case of counting occurences of integers from 1 to 50, tabulate is by far the fastest: about a 150x improvement above <code>table</code>. <a target="_blank" rel="nofollow" href="http://varianceexplained.org/r/board-game-simulation/#fnref:tabulate">↩</a></p>
</li>
<li>
<p>Why did I use reduce from the <code>purrr</code> package rather than the built in <code>Reduce</code>? I like the conciseness of defining each step as <code>~ c(., sum(tail(.)) / 6)</code> rather than <code>function(x, y) c(., sum(tail(.)) / 6)</code>. Note also that while this step builds up a vector incrementally, which is normally a performance hit in R, this calculation takes about a millisecond so it’s not really worth optimizing. <a target="_blank" rel="nofollow" href="http://varianceexplained.org/r/board-game-simulation/#fnref:reduce">↩</a></p>
</li>
</ol>
</div>
<p><div><div><div><div data-width="55" data-show-faces="false" data-send="true" data-layout="box_count" data-height="62" data-href="https://www.r-bloggers.com/the-deadly-board-game-puzzle-efficient-simulation-in-r/"></div></div><div><a data-url="https://www.r-bloggers.com/the-deadly-board-game-puzzle-efficient-simulation-in-r/" data-counturl="https://www.r-bloggers.com/the-deadly-board-game-puzzle-efficient-simulation-in-r/" data-count="vertical" data-text="The ‘deadly board game’ puzzle: efficient simulation in R" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://varianceexplained.org/r/board-game-simulation/"> Variance Explained</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/?s=data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/?s=Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/?s=map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/?s=animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/?s=RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/?s=sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/?s=LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/?s=SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/?s=eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/?s=git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/?s=hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/?s=Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/?s=PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/?s=time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/?s=trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/-mhD_j_jCSg"></a></div>