Return to Revisions

This is related to the Coupon collector's problem as noted in the comments.

Building off of this post, the probability of observing $k$ unique letters in $m$ random samples from an alphabet of size $n$ is:

$\big\{\!{m\!\atop{k}}\big\}\binom{n}{k}\frac{k!}{n^m}=\big\{\!{m\!\atop{k}}\big\}\frac{n!}{n^m(n-k)!}$

Where $\big\{\!{m\!\atop{k}}\big\}$ is the Stirling number of the second kind.

For large $m$, $\ln\big(\big\{\!{m\!\atop{k}}\big\}\big)$ can be approximated.

Here is an R function that returns the probability of every $k$:

library(copula)

coupons1 <- function(n, m) {
  l <- min(m, n)
  k <- 1:l
  if (m < 200) {
    logS <- log(Stirling2.all(n)[k])
  } else {
    # estimate the log Stirling numbers
    v <- n/k
    G <- 1/v
    vexpv <- v/exp(v)
    for (i in 1:5) G <- G - (G - (vexpG <- vexpv*exp(G)))/(1 - vexpG) # Newton's method
    if (l == n) G[l] <- 1
    logS <- (log(v - 1) - log(v*(1 - G)))/2 + (n - k)*(log(v - 1) - log(v - G)) + n*log(k) - k*log(n) + k*(1 - G) + lgamma(n + 1) - lgamma(k + 1) - lgamma(n - k + 1)
  }
  exp(logS + lgamma(n + 1) - lgamma(n - k + 1) - m*log(n))
}

The probability of $k=1,2...43000$ for $n=m=43000$:

system.time(k1 <- coupons1(43e3, 43e3))
#>    user  system elapsed 
#>    0.02    0.00    0.01

plot(26900:27500, k1[26900:27500], xlab = "k", ylab = "p(k)", col = "blue")

Comparing that result to a brute-force approach:

Rcpp::cppFunction("
  NumericVector coupons2(const int& n, const int& m) {
    int maxk;
    int n1 = n - 1;
    if (n > m) {
      maxk = m;
    } else {
      maxk = n;
    }
    
    NumericVector k (maxk);
    k(0) = 1;
    
    for (int i = 1; i < maxk; i++) {
      for (int j = i; j > 0; j--) {
        k(j) = (k(j - 1)*(n - j) + k(j)*(j + 1))/n;
      }
      k(0) = k(0)/n;
    }
    
    for (int i = maxk; i < m; i++) {
      for (int j = n1; j > 0; j--) {
        k(j) = (k(j - 1)*(n - j) + k(j)*(j + 1))/n;
      }
      k(0) = k(0)/n;
    }
    
    return k;
  }
")

system.time(k2 <- coupons2(43e3, 43e3))
#>    user  system elapsed 
#>   12.29    0.00   12.31

The relative error from using the Stirling number approximation is small for a large $m$.

max(abs(k1[26900:27500] - k2[26900:27500])/k1[26900:27500])
#> [1] 8.288002e-07

points(26900:27500, k2[26900:27500], col = "orange", pch = 20)
legend("topright", legend = c("k1", "k2"), col = c("blue", "orange"), pch = c(1, 20))