While the last riddle on The Riddler 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 solution
![\mathbb{E}[Y]=n(1-\frac{1}{n})^m,](https://s0.wp.com/latex.php?latex=%5Cmathbb%7BE%7D%5BY%5D%3Dn%281-%5Cfrac%7B1%7D%7Bn%7D%29%5Em%2C&bg=000000&fg=B0B0B0&s=0&c=20201002)
[which still has a link with e in that the fraction of empty bins converges to e⁻¹ when n=m], this led me to some more involved investigation on the distribution of Y. While it can be shown directly that the probability that k bins are non-empty is
with an R representation by
[sourcecode language=”r” gutter=”false”]
miss<-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)}
[/sourcecode]
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 I came upon this formula
![\mathbb{E}[Y^k]=\sum_{i=1}^k {k \choose i} i! S(k,i) \left( 1-\frac{i}{n}\right)^m](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&fg=B0B0B0&s=0&c=20201002)
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:
[sourcecode language=”r” gutter=”false”]
diss<-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))}
[/sourcecode]
that I still checked by raw simulations from the multinomial
[sourcecode language=”r” gutter=”false”]
zample<-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)}
[/sourcecode]
![\mathbb{E}[|X|]=\mu(2\text{P}_{\mu,\nu}(X>0)-1)\qquad\qquad\\ \qquad\qquad+2 f(0,\nu)\dfrac{\nu}{\nu-1}\left[1+\mu^2/\nu\right]^{-(\nu-1)/2}](https://s0.wp.com/latex.php?latex=%5Cmathbb%7BE%7D%5B%7CX%7C%5D%3D%5Cmu%282%5Ctext%7BP%7D_%7B%5Cmu%2C%5Cnu%7D%28X%3E0%29-1%29%5Cqquad%5Cqquad%5C%5C+%5Cqquad%5Cqquad%2B2+f%280%2C%5Cnu%29%5Cdfrac%7B%5Cnu%7D%7B%5Cnu-1%7D%5Cleft%5B1%2B%5Cmu%5E2%2F%5Cnu%5Cright%5D%5E%7B-%28%5Cnu-1%29%2F2%7D&bg=000000&fg=B0B0B0&s=0&c=20201002)
is the density of the standard Student’s distribution. (I also checked by simulation. The derivation is 
