Oszacuj wielkość populacji, z której pobierana jest próbka, na podstawie liczby powtórzeń obserwacji

Powiedzmy, że mam populację 50 milionów unikalnych rzeczy i pobieram 10 milionów próbek (z wymianą) ... Pierwszy wykres, do którego załączyłem pokazuje, ile razy próbkuję tę samą „rzecz”, co jest stosunkowo rzadkie populacja jest większa niż moja próba.

Jeśli jednak moja populacja liczy tylko 10 milionów rzeczy, a ja pobieram 10 milionów próbek, jak pokazuje drugi wykres, częściej próbuję tę samą rzecz wielokrotnie.

Moje pytanie brzmi - z mojej tabeli obserwacji częstotliwości (dane na wykresach słupkowych) można uzyskać oszacowanie pierwotnej wielkości populacji, gdy jest ona nieznana? Byłoby wspaniale, gdybyś mógł wskazać, jak to zrobić w R.

Jak tam Garvan?

Problem polega na tym, że nie wiemy, ile zaobserwowano zliczeń zerowych. Musimy to oszacować. Klasyczną procedurą statystyczną dla takich sytuacji jest algorytm Expectation-Maximization.

Prosty przykład:

Załóżmy, że czerpiemy z nieznanej populacji (1 000 000) ze stałą poissona równą 0,2.

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

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

Ale nie obserwujemy zera. Zamiast tego obserwujemy to:

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

table

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

Możliwe zaobserwowane częstotliwości

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

Zainicjuj średnią rozkładu Poissona - po prostu zgadnij (wiemy, że tutaj jest 0.2).

lambda <- 1 

1. Oczekiwanie - rozkład Poissona

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

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

3. Druga iteracja

   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

Teraz iteruj aż do konwergencji:

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

Nasza szacunkowa liczba ludności wynosi 1003774, a wskaźnik Poissona szacowany jest na 0,1994473 - jest to szacunkowy odsetek populacji, z której pobrano próbki. Głównym problemem, jaki będziesz mieć w typowych problemach biologicznych, z którymi masz do czynienia, jest założenie, że współczynnik Poissona jest stały.

Przepraszamy za długi post - ta wiki nie nadaje się do kodu R.

To brzmi jak forma „mark and recapture”, czyli „capture-recapture”, dobrze znana technika w ekologii (i niektórych innych dziedzinach, takich jak epidemiologia). Nie moja dziedzina, ale artykuł Wikipedii na temat oznaczania i przechwytywania wygląda rozsądnie, chociaż twoja sytuacja nie jest tą, do której ma zastosowanie wyjaśniona tam metoda Lincolna-Petersena.

Myślę, że shabbychef jest właściwą ścieżką dla twojej sytuacji, ale użycie rozkładu Poissona do przybliżenia dwumianu prawdopodobnie uprościłoby sprawę i powinno być bardzo dobrym przybliżeniem, jeśli populacja jest bardzo duża, jak w twoich przykładach. Wydaje mi się, że uzyskanie wyraźnego wyrażenia dla oszacowania maksymalnego prawdopodobieństwa wielkości populacji powinno być dość proste (patrz ponownie np. Wikipedia ), chociaż nie mam teraz czasu, aby dopracować szczegóły.

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