Randomly Sample M samples from N numbers with replacement, how to estimate N? Can you estimate $N$ with MLE or method of moment or whatever strategy?

*

*$N$ numbered balls are in a bag. $N$ is unknown.

*Pick a ball uniformly at random, record its number, replace it, shuffle.

*After $M$ samples, of which we noticed $R$ repeated numbers, how can we estimate the value of $N$?


Attempt:
If there are $n$ elements of the set then the probability that $x$ have been selected after a sample of $m$ (with replacement) is
$$\frac{S_2(x,m) \; n!}{n^m \; (n-x)!} $$
And I got stuck. Any idea?
 A: This is a standard statistical inference problem involving the classical occupancy distribution (see e.g., O'Neill 2019).  Since $R$ is the number of repeated balls, the number of distinct balls selected in the sample is given by:
$$K = N-R \ \sim \ \text{Occ}(N, M).$$
The probability mass function for this random variable is:
$$p(K=k|N,M) = \frac{(N)_k \cdot S(M,k)}{N^M} \cdot \mathbb{I}(1 \leqslant k \leqslant \min(M,N)),$$
where the values $S(M,k)$ are the Stirling numbers of the second kind and $(N)_k$ are the falling factorials.  The classical occupancy distribution has been subject to a great deal of analysis in the statistical literature, including analysis of statistical inference for the size parameter $N$ (see e.g., Harris 1968).  The form of this distribution and its moments is known, so deriving the MLE or MOM estimators is a relatively simple task.

Maximum-likelihood estimator (MLE): Since the size parameter is an integer, we can find the MLE using discrete calculus.  For any value $1 \leqslant k \leqslant \min(M,N)$ the forward difference of the probability mass function with respect to $N$ can be written as:
$$\begin{align}
\Delta_N p(k) 
&\equiv p(K=k|N+1,M) - p(K=k|N,M) \\[10pt]
&= \frac{(N+1)_k \cdot S(M,k)}{(N+1)^M} - \frac{(N)_k \cdot S(M,k)}{N^M} \\[6pt]
&= S(M,k) \bigg[ \frac{(N+1)_k}{(N+1)^M} - \frac{(N)_k}{N^M} \bigg] \\[6pt]
&= S(M,k) \cdot \frac{(N)_{k}}{(N+1)^M} \bigg[ \frac{N+1}{N-k+1} - \Big( \frac{N+1}{N} \Big)^M \ \bigg] \\[6pt]
\end{align}$$
Thus, if we observe $K=k$ then the maximum-likelihood-estimator (MLE) is given by:
$$\hat{N}_\text{MLE} = \max \bigg \{ N \in \mathbb{N} \ \Bigg| \ \frac{N+1}{N-k+1} < \Big( \frac{N+1}{N} \Big)^M \bigg \}.$$
(There may be cases where the MLE is not unique, since we can also use the $\leqslant$ instead of $<$ in the inequality in this equation.)  Here is a simple function in R to compute the MLE and an example when the input values are fairly large.
MLE.Occ.n <- function(m, k) {
  n <- k
  while ((n+1)/(n-k+1) >= (1+1/n)^m) { n <- n+1 }
  n }

MLE.Occ.n(m = 1000, k = 649)
[1] 1066


Estimation using method-of-moments: The first four moments of the classical occupancy distribution are given in O'Neill (2019) (Section 2).  The expected number of different balls is:
$$\mathbb{E}(K) = N \Bigg[ 1 - \Big( 1-\frac{1}{N} \Big)^M \Bigg].$$
Thus, if we observe $K=k$ then the method-of-moments estimator will approximately solve the implicit equation:
$$\log \hat{N}_\text{MOM}^* - \log k + \text{log1mexp} \Bigg[ - M \log \Big( 1-\frac{1}{\hat{N}_\text{MOM}^*} \Big) \Bigg] = 0.$$
You can solve this equation numerically to obtain a real value $\hat{N}_\text{MOM}^*$ and then use one of the two surrounding integers as $\hat{N}_\text{MOM}$ (these each give slight over- and under-estimates for the true expected value and you can then pick between these using some appropriate method --- e.g., rounding to the nearest integer).  Here is a function in R to compute the method-of-moment estimator.  As can be seen, it gives the same result as the MLE in the present example.
MOM.Occ.n <- function(m, k) {
  FF     <- function(n) { log(n) - log(k) + VGAM::log1mexp(-m*log(1-1/n)) }
  UPPER  <- m*k/(m-k)
  n.real <- uniroot(f = FF, lower = k, upper = UPPER)$root
  round(n.real, 0) }

MOM.Occ.n(m = 1000, k = 649)
[1] 1066

A: I think your likelihood expression has reversed $x=R$ and $m=M$ in $S_2(x,m)$ but no matter - this is a constant with respect to $N$ and so can be ignored.  What you want is the integer $N$ which maximises $\frac{N!}{N^M \; (N-R)!}$.  So you want the largest $N$ where $\frac{N!}{N^M \; (N-R)!} \ge \frac{(N-1)!}{(N-1)^M \; (N-1-R)!} $, i.e. where $N\left(\frac{N-1}{N}\right)^M\ge N-R$, though I doubt this has a simple closed form for $N$.
Another possible approach using a method of moments might be to consider a particular ball so the probability it is never selected is $\left(\frac{N-1}{N}\right)^M$, and the expected number of balls never selected is $N\left(\frac{N-1}{N}\right)^M$ and the expected number selected at least once is $N - N\left(\frac{N-1}{N}\right)^M$, If you see $R$ distinct balls from $M$ attempts then you could try to solve $R= N - N\left(\frac{N-1}{N}\right)^M$ for $N$.  This is essentially the same equation as the likelihood approach, though without the rounding down.
Solving this would not be easy, but in some cases you could use the approximation $\left(\frac{N-1}{N}\right)^M \approx e^{-M/N}$ in which case you might consider $$\hat N\approx \dfrac{M}{\frac{M}{R}+ W\left(-\frac MRe^{-M/R}\right)}$$ where $W$ is the Lambert W function. (When $M \gg R$ the denominator is almost $\frac MR$ so $\hat N$ is very slightly more than $R$, as one might expect.)
As an illustration, if $M=100$ and $R=50$ then direct calculation would eventually give you $\hat N \approx 62.41$ while the suggested approximation could give you $\hat N\approx 62.75$.  The likelihood approach would say $\hat N \le 62.41$ so round this down to $\hat N =62$.
A: I think you would need another constraint.  As described, it would only be possible to estimate a lower bound on the number.   There could be any number of balls.
I think you needed to specify that each ball in the bag has a unique number.
