<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="http://www.r-bloggers.com/occupancy-rules/"></div></div><div><a data-url="http://www.r-bloggers.com/occupancy-rules/" data-counturl="http://www.r-bloggers.com/occupancy-rules/" data-count="horizontal" data-text="occupancy rules" 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="https://xianblog.wordpress.com/2016/05/23/occupancy-rules/"> R – Xi'an's Og</a></strong>, and kindly contributed to <a rel="nofollow" href="http://www.r-bloggers.com/">R-bloggers)</a> </div>
<p><strong>W</strong>hile the last riddle on <a target="_blank" href="http://fivethirtyeight.com/features/can-you-solve-this-elevator-button-puzzle/">The Riddler</a> was rather anticlimactic, namely to find the mean of the number Y of empty bins in a uniform multinomial with n bins and m draws, with <a target="_blank" href="http://math.stackexchange.com/a/1020205/14893">solution</a></p>
<p><img title="\mathbb{E}[Y]=n(1-\frac{1}{n})^m," alt="\mathbb{E}[Y]=n(1-\frac{1}{n})^m," class="latex" src="https://s0.wp.com/latex.php?latex=%5Cmathbb%7BE%7D%5BY%5D%3Dn%281-%5Cfrac%7B1%7D%7Bn%7D%29%5Em%2C&bg=000000&%23038;fg=B0B0B0&%23038;s=0"></p>
<p>[which still has a <a target="_blank" href="https://xianblog.wordpress.com/2016/02/12/the-answer-is-e-what-was-the-question/">link with e</a> in that the fraction of empty bins converges to 1-e⁻¹ when n=m], this led me to some more involved investigation on the distribution of Y. While <a target="_blank" href="http://math.stackexchange.com/a/404055/14893">it can be shown</a> directly that the probability that k bins are non-empty is</p>
<p><img title="{n \choose k}\sum_{i=1}^k (-1)^{k-i}{k \choose i}(i/n)^m" alt="{n \choose k}\sum_{i=1}^k (-1)^{k-i}{k \choose i}(i/n)^m" class="latex" src="https://s0.wp.com/latex.php?latex=%7Bn+%5Cchoose+k%7D%5Csum_%7Bi%3D1%7D%5Ek+%28-1%29%5E%7Bk-i%7D%7Bk+%5Cchoose+i%7D%28i%2Fn%29%5Em&bg=000000&%23038;fg=B0B0B0&%23038;s=0"></p>
<p>with an R representation by</p>
<pre>
miss&lt;-function(n,m){
p=rep(0,n)
for (k in 1:n)
 p[k]=choose(n,k)*sum((-1)^((k-1):0)*choose(k,1:k)*(1:k)^m)
return(rev(p)/n^m)}
</pre>
<p>I wanted to take advantage of the moments of Y, since it writes as a sum of n indicators, counting the number of empty cells. However, the higher moments of Y are not as straightforward as its expectation and I struggled with the representation until <a target="_blank" href="http://mathoverflow.net/a/108665/64496">I came upon this formula</a></p>
<p><img title="\mathbb{E}[Y^k]=\sum_{i=1}^k {k \choose i} i! S(k,i) \left( 1-\frac{i}{n}\right)^m" alt="\mathbb{E}[Y^k]=\sum_{i=1}^k {k \choose i} i! S(k,i) \left( 1-\frac{i}{n}\right)^m" class="latex" src="https://s0.wp.com/latex.php?latex=%5Cmathbb%7BE%7D%5BY%5Ek%5D%3D%5Csum_%7Bi%3D1%7D%5Ek+%7Bk+%5Cchoose+i%7D+i%21+S%28k%2Ci%29+%5Cleft%28+1-%5Cfrac%7Bi%7D%7Bn%7D%5Cright%29%5Em&bg=000000&%23038;fg=B0B0B0&%23038;s=0"></p>
<p>where S(k,i) denotes the Stirling number of the second kind… Or i!S(n,i) is the number of surjections from a set of size n to a set of size i. Which leads to the distribution of Y by inverting the moment equations, as in the following R code:</p>
<pre>
diss&lt;-function(n,m){
  A=matrix(0,n,n)
  mome=rep(0,n)
  A[n,]=rep(1,n)
  mome[n]=1
  for (k in 1:(n-1)){
   A[k,]=(0:(n-1))^k
   for (i in 1:k)
     mome[k]=mome[k]+factorial(i)*as.integer(Stirling2(n,i))*
     (1-(i+1)/n)^m*factorial(k)/factorial(k-i-1)}
  return(solve(A,mome))}
</pre>
<p>that I still checked by raw simulations from the multinomial</p>
<pre>
zample&lt;-function(n,m,T=1e4){
  x=matrix(sample(1:n,m*T,rep=TRUE),nrow=T)
  x=sapply(apply(x,1,unique),length)
  return(n-x)}
</pre>
<p>Filed under: <a target="_blank" href="https://xianblog.wordpress.com/category/kids/">Kids</a>, <a target="_blank" href="https://xianblog.wordpress.com/category/statistics/r-statistics/">R</a>, <a target="_blank" href="https://xianblog.wordpress.com/category/statistics/">Statistics</a> Tagged: <a target="_blank" href="https://xianblog.wordpress.com/tag/moment-derivation/">moment derivation</a>, <a target="_blank" href="https://xianblog.wordpress.com/tag/moments/">moments</a>, <a target="_blank" href="https://xianblog.wordpress.com/tag/multinomial-distribution/">multinomial distribution</a>, <a target="_blank" href="https://xianblog.wordpress.com/tag/occupancy/">occupancy</a>, <a target="_blank" href="https://xianblog.wordpress.com/tag/r/">R</a>, <a target="_blank" href="https://xianblog.wordpress.com/tag/stack-exchange/">Stack Exchange</a>, <a target="_blank" href="https://xianblog.wordpress.com/tag/stirling-number/">Stirling number</a>, <a target="_blank" href="https://xianblog.wordpress.com/tag/surjection/">surjection</a> <a target="_blank" rel="nofollow" href="http://feeds.wordpress.com/1.0/gocomments/xianblog.wordpress.com/30620/"><img alt="" src="http://feeds.wordpress.com/1.0/comments/xianblog.wordpress.com/30620/" border="0"></a> <img height="1" alt="" width="1" src="https://pixel.wp.com/b.gif?host=xianblog.wordpress.com&blog=5051449&%23038;post=30620&%23038;subd=xianblog&%23038;ref=&%23038;feed=1" border="0"></p>
<p><div><div><div><div data-width="55" data-show-faces="false" data-send="true" data-layout="box_count" data-height="62" data-href="http://www.r-bloggers.com/occupancy-rules/"></div></div><div><a data-url="http://www.r-bloggers.com/occupancy-rules/" data-counturl="http://www.r-bloggers.com/occupancy-rules/" data-count="vertical" data-text="occupancy rules" 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="https://xianblog.wordpress.com/2016/05/23/occupancy-rules/"> R – Xi'an's Og</a></strong>.</div>
<hr>
<a rel="nofollow" href="http://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="http://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="http://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="http://www.r-bloggers.com/?s=data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="http://www.r-bloggers.com/?s=Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="http://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="http://www.r-bloggers.com/?s=map">maps</a>, <a title="Animation in R" rel="nofollow" href="http://www.r-bloggers.com/?s=animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="http://www.r-bloggers.com/?s=RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="http://www.r-bloggers.com/?s=sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="http://www.r-bloggers.com/?s=LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="http://www.r-bloggers.com/?s=SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="http://www.r-bloggers.com/?s=eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="http://www.r-bloggers.com/?s=git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="http://www.r-bloggers.com/?s=hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="http://www.r-bloggers.com/?s=Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="http://www.r-bloggers.com/?s=PCA">PCA</a>, <a title="Time series" rel="nofollow" href="http://www.r-bloggers.com/?s=time+series">time series</a>, <a title="finance trading" rel="nofollow" href="http://www.r-bloggers.com/?s=trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/yRvHtKyFoFI"></a></div>