How to estimate total number of different results for a stochastic event?

Question

Let $X$ be a random variable that has $n$ possible values $x_1,x_2,\dots,x_n$, and that $P(X=x_i)=\frac 1n ,\forall x=1\dots n$.

Note that $n$ and $x_i,i=1\dots n$ are all unknown and $x_i$ is unordered, but with any two results $a_i, a_j\in\{x_1\dots x_n\}$, you can know whether they are the same value.

Assume the random event has been observed $m$ times, and $p$ different values $a_1, a_2, \dots, a_p$ have been observed, with each value appearing $c_1, c_2, \dots, c_p$ times. Apparently there are $c_i>0$ and $\sum c_i=m$. Assume $p$ is significantly smaller than $n$ ($p<<n$). Given the knowledge that each $x_i$ have an equal chance of $\frac 1n$ of appearing for any single test (observation) and the array $\{c_i\}$, how to make an unbiased estimation of $n$?

---

A version that's easier to understand:

You're playing a slot machine with only one slot. You don't know how many different items are there on the slot, but you know that the machine is "uncheated" and each possible item has an equal chance of appearing on every roll. You've rolled the slot machine for $m$ times and seen $p$ different items. You count your results and see that each item has appeared for $c_i$ times. You know that the slot machine has significantly more items available than you have seen. How do you make an unbiased estimate of the total number of items on the slot roll?

---

I'm aware of the German tank problem but found mine vastly different in many ways, such as the outcomes of the random event being unordered and the statistics of observed results matters, so I couldn't apply the German tank model and work this out by myself.

I'm also assuming that the pattern is "typical", such that very few items appear for some more times (above 3 or 4), and slightly more items appearing twice or thrice, and the vast majority appearing only once.

Disclaimer: This is not in any form related to schoolwork. I came up with this question while playing computer games.

Answer 1

Suppose 4 events are observed with frequency vector $c=(1,1,2)$ as in @quasi 's example and that there are really 5 unique objects labeled A, B, C, D, and E. The probability of observing one A, two D's, and one E is given by the multinomial probability mass function:

$$\text{multinomial} = \frac{4! \left(\frac{1}{5}\right)^4}{1! 0! 0! 2! 1!}=\frac{12}{625}$$

But because we don't know if we've seen A, D, and E or B, C, and D or A, B, and E, etc. we need to multiply that probability by the number of possible arrangements of the selected objects. To do so we look at the frequency of the frequencies. We have the "true" frequencies of (1, 1, 2, 0, 0). There are 2 objects with frequency 1 and 1 object with frequency 2, and 2 objects with frequency 0. That frequency of frequency vector is $f = (2,1,2)$. The possible number of arrangements is

$$\text{multiplier} = \frac{5!}{2! 1! 2!}=30$$

So the probability of the observed frequencies $c=(1,1,2)$ is multinomial*multiplier = (12/625)*30 = 72/125 = 0.576.

You go through this process for $n = 3, 4, 5, 6, \ldots$ and choose the value of $n$ that maximizes the probability of the observed frequencies.

Some Mathematica code to do this for a proposed set of observed frequencies follows:

prob[c_] := (Total[c]!/((c!) /. List -> Times)) (1/Length[c])^Total[c] * 
  (Length[c]!/((Tally[c][[All, 2]]!) /. List -> Times))

{3, prob[{1, 1, 2}] // N}
(* {3, 0.444444} *)
{4, prob[{1, 1, 2, 0}] // N}
(* {4, 0.5625} *)
{5, prob[{1, 1, 2, 0, 0}] // N}
(* {5, 0.576} *)
{6, prob[{1, 1, 2, 0, 0, 0}] // N}
(* {6, 0.555556} *)
{7, prob[{1, 1, 2, 0, 0, 0, 0}] // N}
(* {7, 0.524781} *)
{8, prob[{1, 1, 2, 0, 0, 0, 0, 0}] // N}
(* {8, 0.492188} *)

We see that $n=5$ maximizes probability of observing $c=(1,1,2)$.

That is the process for determinine the maximum likelihood estimate given a particular set of observed frequencies. What is also important is knowing the distribution of the maximum likelihood estimator given the sample size ($m$) and the number of unique elements in the population ($n$).

Because the maximum likelihood estimate is $\infty$ when all of the observed frequencies are 1, the maximum likelihood estimator has no mean and therefore can't be unbiased (as you mentioned that unbiasedness was important to you). That doesn't mean that there aren't unbiased estimators but just that using maximum likelihood won't achieve that.

Here is some Mathematica code to obtain the distribution of the maximum likelihood estimator of $n$ given the sample size $m$. First, define a few functions to obtain the possible samples, probabilities, and maximum likelihood estimates:

(* List of possible observed frequencies given sample size and number of items in population *)
ss[m_, n_] := 
 If[Length[#] < n, Join[#, ConstantArray[0, n - Length[#]]], #] & /@ IntegerPartitions[m, {1, n}]

(* Probability of observing a particular set of n frequencies *)
prob[c_] := (Total[c]!/((c!) /. List -> Times)) (1/Length[c])^Total[c] * 
  (Length[c]!/((Tally[c][[All, 2]]!) /. List -> Times))

(* Maximum likelihood estimate of n given observed frequency counts *)    
mle[c_] := Module[{n0},
  n0 = Length[c];
  If[Total[c] == Length[c], \[Infinity],
   Sort[Join[{{n0, prob[c] // N}},
      Table[{i, prob[Join[c, ConstantArray[0, i - n0]]] // N}, {i, n0 + 1, 500}]],
     #1[[2]] > #2[[2]] &][[1, 1]]]]

(Note that the mle function only allows a maximum value of $n$ being 500. That maximum can be increased if 500 is ever reached.) Now use the functions to obtain the distribution of the maximum likelihood estimator:

m = 10; (* Sample size *)
n = 20;  (* Number of items in population *)

(* Determine distribution of the maximum likelihood estimator given m and n *)
data = Transpose[{mle[#] & /@ IntegerPartitions[m, {1, n}], 
   prob[#] & /@ ss[m, n]}];
g = GatherBy[data, #[[1]] &];
dist = {#[[1, 1]], Total[#[[All, 2]]] // N} & /@ g;
TableForm[dist, TableHeadings -> {None, {"MLE", "Probability"}}]

The estimation problem you describe is related to capture/recapture statistical procedures and so likely this is a well-known topic (just not well-known to me). A Bayesian approach might be fruitful if you can characterize what you think about the possible values of $n$ as a probability distribution.

Answer 2

Here is an approach based on Example 4.11 in Combinatorial Methods in Discrete Distributions by C. A. Charalambides (2005), Wiley Series in Probability and Statistics.

When $m > n$ (sample size greater than the number of unique individuals in the population) there exists a uniformly minimum-variance unbiased estimator of $n$. For $m < n$ the estimator is biased. Using the functions in my other answer one can determine the bias and variance. That estimator is

$$\hat{n}=p + \frac{\mathcal{S}_m^{(p-1)}}{\mathcal{S}_m^{(p)}}$$

where $\mathcal{S}_m^{(p)}$ is the Stirling number of the second kind.

An estimate of the variance for a particular observed $p$ and $m$ is given by

$$\hat{v}=\left(\frac{\mathcal{S}_m^{(p-1)}}{\mathcal{S}_m^{(p)}}\right){}^2+\frac{p \mathcal{S}_m^{(p-1)}}{\mathcal{S}_m^{(p)}}-\frac{\mathcal{S}_m^{(p-2)}}{\mathcal{S}_m^{(p)}}$$

For example, if one takes a sample of size 12 ($m=12$) and observes $(c_1,c_2,c_3,c_4,c_5,c_6,c_7,c_8)=(1, 1, 1, 1, 1, 2, 2, 3)$, then $p=8$ and the estimate of $n$ is

$$8+\frac{\mathcal{S}_{12}^{(7)}}{\mathcal{S}_{12}^{(8)}}=\frac{57564}{4819}\approx11.9452$$