How to guess the size of a set? Assume we have a set of unique words and draw a number $n$ of them using simple-random-sampling without replacement independently in each round. We have several rounds and try to guess the set size after each round. The words are drawn from the same set in each round, but are drawn without replacement in each round (all words are distinct in each round). How could this be achieved and how big would the error be?
For example if the set size is smaller than $n$, then we would only need one round.  If the set size is close to $n$, in the second round the number of words that were drawn in the first round would be close to $n$.  So I know that we get some information, but I don't understand how we could efficiently guess the set size.

Using the following notation:
$s_i$: the number of samples for the $i^{th}$ round
$k_i$: the number of words sampled in the $i^{th}$ round that had not been previously sampled
$m_i=\sum_{j=1}^ik_j$
$\pi_1(n)$: the prior distribution of $n$
We'll also assume $n\geq{s}$; otherwise we can determine $n$ after the first round of sampling. We can update $\pi(n)$ beginning with the second round of sampling:
$\pi_i(n)\propto\frac{\binom{m_{i-1}}{s_i-k_i}\binom{n-m_{i-1}}{k_i}}{\binom{n}{s_i}}\pi_{i-1}(n)$
Implemented as an R function:
pn <- function(s, k, prior) {
  l <- length(s)
  m <- cumsum(k[-l])
  s <- s[-1]
  k <- k[-1]
  k1 <- k + 1
  m1 <- m + 1
  s1 <- s + 1
  function(n) pmax(prior(n)*exp(sum(lgamma(c(m1, s1))) - sum(lgamma(c(s1 - k, m1 - s + k, k1))) + colSums(lgamma(outer(-m, n + 1, "+")) + lgamma(outer(-s, n + 1, "+")) - lgamma(outer(-m - k, n + 1, "+"))) - (l - 1)*lgamma(n + 1)), 0, na.rm = TRUE)
}

For example, say $\pi_1(n)\sim{U(8,30)}$, $k=5,3,1$, and $s=5,5,5$.
k <- c(5, 3, 1)
s <- rep(5, 3)

post <- pn(s, k, function(n) 1)
like <- post(8:30)
plot(8:30, like/sum(like), xlab = "n", ylab = "p(n)", col = "blue", pch = 3)


We can verify the results with simulation:
library(parallel)
set.seed(724526144)
nreps <- tabulate(sample(23, 23e6, TRUE), 23)
clust <- makeCluster(detectCores() - 1)
clusterExport(clust, "nreps")
sim <- unlist(parLapply(clust, 1:23, function(i) sum(replicate(nreps[i], all(cumsum(!duplicated(c(replicate(3, sample(i + 7, 5)))))[c(5, 10, 15)] == c(5L, 8L, 9L))))))
stopCluster(clust)
points(8:30, sim/sum(sim), col = "orange")
legend("topright", legend = c("Posterior likelihood", "Simulation"), col = c("blue", "orange"), pch = c(3, 1))


If we had observed $k=5,4,2$ instead of $k=5,3,1$, we can see that the posterior distribution would shift toward larger values of $n$.
k <- c(5, 4, 2)

post <- pn(s, k, function(n) 1)
like <- post(8:30)
plot(8:30, like/sum(like), xlab = "n", ylab = "p(n)", col = "blue", pch = 3)


Again verifying the results with simulation:
set.seed(353678169)
nreps <- tabulate(sample(23, 23e6, TRUE), 23)
clust <- makeCluster(detectCores() - 1)
clusterExport(clust, "nreps")
sim <- unlist(parLapply(clust, 1:23, function(i) sum(replicate(nreps[i], all(cumsum(!duplicated(c(replicate(3, sample(i + 7, 5)))))[c(5, 10, 15)] == c(5L, 9L, 11L))))))
stopCluster(clust)
points(8:30, sim/sum(sim), col = "orange")
legend("topright", legend = c("Posterior likelihood", "Simulation"), col = c("blue", "orange"), pch = c(3, 1))


