Schätzen Sie die Größe einer Population, die beprobt wird, anhand der Anzahl der wiederholten Beobachtungen

Angenommen, ich habe eine Bevölkerung von 50 Millionen Unikaten, und ich nehme 10 Millionen Proben (mit Ersatz) ... Die erste Grafik, die ich beigefügt habe, zeigt, wie oft ich das gleiche "Ding" probiere, was relativ selten ist wie die Bevölkerung ist größer als meine Stichprobe.

Wenn meine Bevölkerung jedoch nur aus 10 Millionen Dingen besteht und ich 10 Millionen Proben nehme, zeige die zweite Grafik, dass ich die gleiche Sache öfter und mehrmals probiere.

Meine Frage ist: Ist es aus meiner Häufigkeitstabelle (den Daten in den Balkendiagrammen) möglich, eine Schätzung der ursprünglichen Populationsgröße zu erhalten, wenn diese unbekannt ist? Und es wäre großartig, wenn Sie einen Hinweis geben könnten, wie dies in R geschehen soll.

Wie geht es dem Garvan?

Das Problem ist, dass wir nicht wissen, wie viele Nullzählungen beobachtet werden. Das müssen wir abschätzen. Ein klassisches statistisches Verfahren für Situationen wie diese ist der Expectation-Maximization-Algorithmus.

Ein einfaches Beispiel:

Angenommen, wir ziehen aus einer unbekannten Population (von 1.000.000) mit einer Poisson-Konstante von 0,2.

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

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

Aber wir beobachten die Nullzählungen nicht. Stattdessen beobachten wir Folgendes:

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

table

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

Mögliche Frequenzen beobachtet

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

Initialisiere den Mittelwert der Poisson-Verteilung - nimm einfach eine Schätzung (wir wissen, dass es hier 0,2 ist).

lambda <- 1 

1. Erwartung - Poissonverteilung

   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. Maximierung

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

3. Zweite Iteration

   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

Nun iteriere bis zur Konvergenz:

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

Unsere Bevölkerungsschätzung beträgt 1003774 und unsere Giftrate wird auf 0,1994473 geschätzt - dies ist der geschätzte Anteil der beprobten Bevölkerung. Das Hauptproblem bei den typischen biologischen Problemen ist die Annahme, dass die Giftrate eine Konstante ist.

Entschuldigung für den langwierigen Beitrag - dieses Wiki ist nicht wirklich für R-Code geeignet.

Dies klingt nach einer Form von „Mark and Recapture“ oder „Capture-Recapture“, einer bekannten Technik in der Ökologie (und einigen anderen Bereichen wie der Epidemiologie). Nicht mein Gebiet, sondern der Wikipedia-Artikel über Mark und Recapture sieht vernünftig aus, obwohl Ihre Situation nicht die ist, auf die die dort erläuterte Lincoln-Petersen-Methode zutrifft.

Ich denke, shabbychef ist der richtige Weg für Ihre Situation, aber die Poisson-Verteilung zur Approximation des Binomials würde die Sache wahrscheinlich etwas einfacher machen und sollte eine sehr gute Approximation sein, wenn die Bevölkerungszahl sehr groß ist, wie in Ihren Beispielen. Ich denke, ein expliziter Ausdruck für die maximale Wahrscheinlichkeitsschätzung der Populationsgröße sollte dann ziemlich einfach sein (siehe z. B. Wikipedia ), obwohl ich momentan keine Zeit habe, die Details zu erarbeiten.

$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$