<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/flea-circus/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://xianblog.wordpress.com/2016/12/08/flea-circus/"> R – Xi'an's Og</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p><a target="_blank" rel="nofollow" href="https://i1.wp.com/xianblog.files.wordpress.com/2011/10/grib1.gif?ssl=1"><img alt="grib" class="wp-image-12730" src="https://i1.wp.com/xianblog.files.wordpress.com/2011/10/grib1.gif?resize=450%2C279&ssl=1" data-recalc-dims="1"></a><strong>A</strong>n old riddle found on <a target="_blank" rel="nofollow" href="http://stats.stackexchange.com/q/2427/7224">X validated</a> asking for Monte Carlo resolution  but originally given on <a target="_blank" rel="nofollow" href="https://projecteuler.net/problem=213">Project Euler</a>:</p>
<blockquote><p><span><em>A 30×30 grid of squares contains 30² fleas, initially one flea per square. When a bell is rung, each flea jumps to an adjacent square at random. What is the expected number of unoccupied squares after 50 bell rings, up to six decimal places?</em></span></p>
</blockquote>
<p>The debate on <a target="_blank" rel="nofollow" href="http://stats.stackexchange.com/q/2427/7224">X validated</a> is whether or not a Monte Carlo resolution is feasible. Up to six decimals, certainly not. But with some lower precision, certainly. Here is a rather basic R code where the 50 steps are operated on the 900 squares, rather than the 900 fleas. This saves some time by avoiding empty squares.</p>
<pre>xprmt=function(n=10,T=50){

 mean=0
 for (t in 1:n){

   board=rep(1,900)
   for (v in 1:T){

    beard=rep(0,900)
    if (board[1]>0){
      poz=c(0,1,0,30)
      ne=rmultinom(1,board[1],prob=(poz!=0))
      beard[1+poz]=beard[1+poz]+ne}
    #
    for (i in (2:899)[board[-1][-899]>0]){
     u=(i-1)%%30+1;v=(i-1)%/%30+1
     poz=c(-(u>1),(u&lt;30),-30*(v>1),30*(v&lt;30))      
     ne=rmultinom(1,board[i],prob=(poz!=0))      
     beard[i+poz]=beard[i+poz]+ne} 
    #     
    if (board[900]>0){
     poz=c(-1,0,-30,0)
     ne=rmultinom(1,board[900],prob=(poz!=0))
     beard[900+poz]=beard[900+poz]+ne}
     board=beard}
   mean=mean+sum(board==0)}
 return(mean/n)}
</pre>
<p>The function returns an empirical average over n replications. With a presumably awkward approach to the borderline squares, since it involves adding zeros to keep the structure the same… Nonetheless, it produces an approximation that is rather close to the approximate expected value, in about 3mn on my laptop.</p>
<pre>> exprmt(n=1e3)
[1] 331.082
> 900/exp(1)
[1] 331.0915
</pre>
<p>Further gains follow from considering only half of the squares, as there are two independent processes acting in parallel. I looked at an alternative and much faster approach using the stationary distribution, with the stationary being the Multinomial (450,(2/1740,3/1740…,4/1740,…,2/1740)) with probabilities proportional to 2 in the corner, 3 on the sides, and 4 in the inside. (The process, strictly speaking, has <em>no stationary distribution</em>, since it is periodic. But one can consider instead the subprocess indexed by even times.) This seems to be the case, though, when looking at the occupancy frequencies, after defining the stationary as:</p>
<pre>inva=function(B=30){
return(c(2,rep(3,B-2),2,rep(c(3,
  rep(4,B-2),3),B-2),2,rep(3,B-2),2))}
</pre>
<p>namely</p>
<pre>> mn=0;n=1e8 #14 clock hours!
> proz=rep(c(rep(c(0,1),15),rep(c(1,0),15)),15)*inva(30)
> for (t in 1:n)
+ mn=mn+table(rmultinom(1,450,prob=rep(1,450)))[1:4]
> mn=mn/n
> mn[1]=mn[1]-450
> mn
     0      1      2     3
166.11 164.92  82.56 27.71
> exprmt(n=1e6) #55 clock hours!!
[1] 165.36 165.69 82.92 27.57</pre>
<p>my original confusion being that the Poisson approximation had not yet taken over… (Of course, computing the first frequency for the stationary distribution does not require any simulation, since it is the sum of the complement probabilities to the power 450, i.e., 166.1069.)</p>
<p>Filed under: <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/category/books/">Books</a>, <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/category/kids/">Kids</a>, <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/category/pictures/">pictures</a>, <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/category/statistics/r-statistics/">R</a>, <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/category/statistics/">Statistics</a> Tagged: <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/tag/coupling/">coupling</a>, <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/tag/cross-validated/">cross validated</a>, <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/tag/fleas/">fleas</a>, <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/tag/leonhard-euler/">Leonhard Euler</a>, <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/tag/markov-chains/">Markov chains</a>, <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/tag/markov-random-field/">Markov random field</a>, <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/tag/monte-carlo-integration/">Monte Carlo integration</a>, <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/tag/occupancy/">occupancy</a>, <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/tag/poisson-approximation/">Poisson approximation</a>, <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/tag/r/">R</a>, <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/tag/random-walk/">random walk</a>, <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/tag/stationarity/">stationarity</a> <a target="_blank" rel="nofollow" href="http://feeds.wordpress.com/1.0/gocomments/xianblog.wordpress.com/32101/"><img alt="" src="http://feeds.wordpress.com/1.0/comments/xianblog.wordpress.com/32101/" border="0"></a> <img alt="" src="https://i1.wp.com/pixel.wp.com/b.gif?resize=1%2C1&ssl=1" border="0" data-recalc-dims="1"></p>
<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://xianblog.wordpress.com/2016/12/08/flea-circus/"> R – Xi'an's Og</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/search/tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/search/data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/search/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/search/ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/search/map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/search/animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/search/sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/search/SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/search/git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/search/hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/search/PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/search/time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/search/trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/-81b3c-re8Q"></a></div>