¿Cómo estimar el número total de resultados diferentes para un evento estocástico?

Question

Dejemos que $X$ sea una variable aleatoria que tiene $n$ valores posibles $x_1,x_2,\dots,x_n$ y que $P(X=x_i)=\frac 1n ,\forall x=1\dots n$ .

Tenga en cuenta que $n$ y $x_i,i=1\dots n$ son todos desconocidos y $x_i$ es desordenado pero con dos resultados cualesquiera $a_i, a_j\in\{x_1\dots x_n\}$ , se puede saber si son del mismo valor.

Supongamos que se ha observado el suceso aleatorio $m$ veces, y $p$ valores diferentes $a_1, a_2, \dots, a_p$ se han observado, apareciendo cada valor $c_1, c_2, \dots, c_p$ tiempos. Aparentemente hay $c_i>0$ y $\sum c_i=m$ . Supongamos que $p$ es significativamente menor que $n$ ( $p<<n$ ). Teniendo en cuenta que cada $x_i$ tienen la misma posibilidad de $\frac 1n$ de aparecer para una sola prueba (observación) y el conjunto $\{c_i\}$ cómo hacer una estimación no sesgada de $n$ ?

---

Una versión más fácil de entender:

Estás jugando a una máquina tragaperras con una sola ranura. No sabe cuántos elementos diferentes hay en la ranura, pero sabe que la máquina no tiene trampa y que cada elemento posible tiene la misma probabilidad de aparecer en cada tirada. Has hecho rodar la máquina tragaperras por $m$ veces y visto $p$ diferentes artículos. Usted cuenta sus resultados y ve que cada artículo ha aparecido por $c_i$ veces. Sabes que la máquina tragaperras tiene muchos más artículos disponibles de los que has visto. ¿Cómo puede hacer una estimación no sesgada del número total de artículos en el rollo de la tragamonedas?

---

Soy consciente de la Problema de los tanques alemanes pero descubrí que el mío es muy diferente en muchos aspectos, como que los resultados del evento aleatorio no están ordenados y que las estadísticas de los resultados observados son importantes, por lo que no pude aplicar el modelo del tanque alemán y resolverlo por mí mismo.

También estoy asumiendo que el patrón es "típico", de manera que muy pocos artículos aparecen algunas veces más (por encima de 3 o 4), y un poco más de artículos que aparecen dos o tres veces, y la gran mayoría que aparecen sólo una vez.

Descargo de responsabilidad: Esto no está relacionado de ninguna manera con el trabajo escolar. Se me ocurrió esta pregunta mientras jugaba con el ordenador.

Answer 1

Supongamos que se observan 4 eventos con el vector de frecuencia $c=(1,1,2)$ como en el ejemplo de @quasi y que realmente hay 5 objetos únicos etiquetados como A, B, C, D y E. La probabilidad de observar un A, dos D y un E viene dada por la función de masa de probabilidad multinomial:

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

Pero como no sabemos si hemos visto A, D y E o B, C y D o A, B y E, etc. tenemos que multiplicar esa probabilidad por el número de disposiciones posibles de los objetos seleccionados. Para ello, miramos la frecuencia de las frecuencias. Tenemos las frecuencias "verdaderas" de (1, 1, 2, 0, 0). Hay 2 objetos con frecuencia 1 y 1 objeto con frecuencia 2, y 2 objetos con frecuencia 0. Ese vector de frecuencias es $f = (2,1,2)$ . El número posible de arreglos es

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

Así que la probabilidad de las frecuencias observadas $c=(1,1,2)$ es multinomio*multiplicador = (12/625)*30 = 72/125 = 0,576.

Usted pasa por este proceso para $n = 3, 4, 5, 6, \ldots$ y elegir el valor de $n$ que maximiza la probabilidad de las frecuencias observadas.

Algunos Mathematica El código para hacer esto para un conjunto propuesto de frecuencias observadas es el siguiente:

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} *)

Vemos que $n=5$ maximiza la probabilidad de observar $c=(1,1,2)$ .

Es el proceso para determinar la estimación de máxima verosimilitud dado un conjunto particular de frecuencias observadas. También es importante conocer la distribución del estimador de máxima verosimilitud dado el tamaño de la muestra ( $m$ ) y el número de elementos únicos en la población ( $n$ ).

Porque la estimación de máxima verosimilitud es $\infty$ cuando todas las frecuencias observadas son 1, el estimador de máxima verosimilitud no tiene media y, por lo tanto, no puede ser insesgado (como usted mencionó que la insesgadez era importante para usted). Eso no significa que no haya estimadores insesgados, sino que el uso de la máxima verosimilitud no lo conseguirá.

Aquí hay algunos Mathematica para obtener la distribución del estimador de máxima verosimilitud de $n$ dado el tamaño de la muestra $m$ . En primer lugar, se definen algunas funciones para obtener las muestras posibles, las probabilidades y las estimaciones de máxima verosimilitud:

(* 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]]]]

(Tenga en cuenta que el mle sólo permite un valor máximo de $n$ siendo 500. Ese máximo puede aumentarse si se alcanza 500). Ahora utilice las funciones para obtener la distribución del estimador de máxima verosimilitud:

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"}}]

El problema de estimación que describes está relacionado con los procedimientos estadísticos de captura/recaptura, por lo que es probable que se trate de un tema bien conocido (sólo que yo no lo conozco). Un enfoque bayesiano podría ser fructífero si puede caracterizar lo que piensa sobre los posibles valores de $n$ como una distribución de probabilidad.

Answer 2

Este es un enfoque basado en el ejemplo 4.11 de Métodos combinatorios en distribuciones discretas de C. A. Charalambides (2005), Wiley Series in Probability and Statistics.

Cuando $m > n$ (tamaño de la muestra mayor que el número de individuos únicos en la población) existe un estimador insesgado de varianza mínima uniforme de $n$ . Para $m < n$ el estimador está sesgado. Utilizando las funciones de mi otra respuesta se puede determinar el sesgo y la varianza. Ese estimador es

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

donde $\mathcal{S}_m^{(p)}$ es el número de Stirling del segundo tipo.

Una estimación de la varianza para un determinado observado $p$ y $m$ viene dada por

$$\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)}}$$

Por ejemplo, si se toma una muestra de tamaño 12 ( $m=12$ ) y observa $(c_1,c_2,c_3,c_4,c_5,c_6,c_7,c_8)=(1, 1, 1, 1, 1, 2, 2, 3)$ entonces $p=8$ y la estimación de $n$ es

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