Estimer la taille d'une population échantillonnée par le nombre d'observations répétées

Disons que j'ai une population de 50 millions de choses uniques, et je prends 10 millions d'échantillons (avec remplacement) ... Le premier graphique que j'ai joint montre combien de fois j'échantillonne la même "chose", ce qui est relativement rare que le la population est plus grande que mon échantillon.

Cependant, si ma population n'est que de 10 millions de choses et que je prends 10 millions d'échantillons, comme le montre le deuxième graphique, j'échantillonnerai plus souvent la même chose plusieurs fois.

Ma question est - à partir de mon tableau de fréquence des observations (les données des graphiques à barres) est-il possible d'obtenir une estimation de la taille de la population d'origine quand elle est inconnue? Et ce serait formidable si vous pouviez fournir un pointeur sur la façon de procéder à ce sujet dans R.

Comment va le Garvan?

Le problème est que nous ne savons pas combien de dénombrements zéro sont observés. Nous devons l'estimer. L'algorithme Expectation-Maximization est une procédure statistique classique pour des situations comme celle-ci.

Un exemple simple:

Supposons que nous puisons dans une population inconnue (de 1 000 000) avec une constante de poisson de 0,2.

counts <- rpois(1000000, 0.2)
table(counts)

     0      1      2      3      4      5
818501 164042  16281   1111     62      3

Mais nous n'observons pas le compte zéro. Au lieu de cela, nous observons ceci:

table <- c("0"=0, table(counts)[2:6])

table

     0      1      2      3      4      5
     0 164042  16281   1111     62      3

Fréquences possibles observées

k <- c("0"=0, "1"=1, "2"=2, "3"=3, "4"=4, "5"=5)

Initialisez la moyenne de la distribution de Poisson - faites une supposition (nous savons que c'est 0,2 ici).

lambda <- 1 

1. Attente - Distribution de Poisson

   P_k <- lambda^k*exp(-lambda)/factorial(k)
   P_k
                 0           1           2           3           4           5
   0.367879441 0.367879441 0.183939721 0.061313240 0.015328310 0.003065662  
   n0 <- sum(table[2:6])/(1 - P_k[1]) - sum(table[2:6])
   
   
   n0
          0
   105628.2     
   table[1] <-  105628.2

2. Maximisation

   lambda_MLE <- (1/sum(table))*(sum(table*k))        
   lambda_MLE        
   [1] 0.697252        
   lambda <- lambda_MLE

3. Deuxième itération

   P_k <- lambda^k*exp(-lambda)/factorial(k)        
   n0 <- sum(table[2:6])/(1 - P_k[1]) - sum(table[2:6])       
   table[1] <-  n0 
   lambda <- (1/sum(table))*(sum(table*k))
   
   
    population lambda_MLE
   
   
   [1,]   361517.1  0.5537774

Maintenant, itérez jusqu'à convergence:

for (i in 1:200) {  
P_k <- lambda^k*exp(-lambda)/factorial(k)  
n0 <- sum(table[2:6])/(1 - P_k[1]) - sum(table[2:6])
table[1] <-  n0
lambda <- (1/sum(table))*(sum(table*k))
}
cbind( population = sum(table), lambda_MLE)
     population lambda_MLE
[1,]    1003774  0.1994473

Notre estimation de la population est de 1003774 et notre taux de poisson est estimé à 0,1994473 - c'est la proportion estimée de la population échantillonnée. Le principal problème que vous rencontrerez dans les problèmes biologiques typiques que vous rencontrez est l'hypothèse que le taux de poisson est constant.

Désolé pour le message de longue haleine - ce wiki n'est pas vraiment adapté au code R.

Cela ressemble à une forme de «marquage et recapture» alias «capture-recapture», une technique bien connue en écologie (et dans certains autres domaines tels que l'épidémiologie). Ce n'est pas mon domaine, mais l'article de Wikipédia sur la marque et la recapture semble raisonnable, bien que votre situation ne soit pas celle à laquelle la méthode Lincoln – Petersen expliquée ici s'applique.

Je pense que shabbychef est l'une des bonnes pistes pour votre situation, mais l'utilisation de la distribution de Poisson pour approximer le binôme simplifierait probablement les choses et devrait être une très bonne approximation si la taille de la population est très grande, comme dans vos exemples. Je pense que l'obtention d'une expression explicite pour l'estimation du maximum de vraisemblance de la taille de la population devrait alors être assez simple (voir par exemple Wikipedia à nouveau ), bien que je n'ai pas le temps de travailler sur les détails pour le moment.

$n$$k$$k$$P = \frac{1}{k}$$m$$m$$n$${n \choose m} P^m (1-P)^{n-m}$$n$$n$$k$$m$$(1-P)$$1$

$P_m$$m$$P_{m} / P_{m+1}$$(k-1)\frac{m+1}{n-m}$$k$