Estimar o tamanho de uma população amostrada pelo número de observações repetidas

13

Digamos que tenho uma população de 50 milhões de coisas únicas e colho 10 milhões de amostras (com substituição) ... O primeiro gráfico anexado mostra quantas vezes eu coleciono a mesma "coisa", o que é relativamente raro como o população é maior que minha amostra.

No entanto, se minha população é de apenas 10 milhões de coisas e eu colho 10 milhões de amostras, como o segundo gráfico mostra, vou experimentar com mais frequência a mesma coisa repetidas vezes.

Minha pergunta é: da minha tabela de observações de frequência (os dados nos gráficos de barras) é possível obter uma estimativa do tamanho da população original quando é desconhecido? E seria ótimo se você pudesse fornecer um ponteiro sobre como fazer isso em R.

texto alternativo

Aaron Statham
fonte

Respostas:

10

Como está o Garvan?

O problema é que não sabemos quantas contagens zero são observadas. Temos que estimar isso. Um procedimento estatístico clássico para situações como essa é o algoritmo Expectation-Maximization.

Um exemplo simples:

Suponha que extraímos de uma população desconhecida (de 1.000.000) com uma 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

Mas não observamos a contagem zero. Em vez disso, observamos o seguinte:

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

table

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

Possíveis frequências observadas

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

Inicialize a média da distribuição de Poisson - basta adivinhar (sabemos que é 0,2 aqui).

lambda <- 1 
  1. Expectativa - Distribuição 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. Maximização

    lambda_MLE <- (1/sum(table))*(sum(table*k))        
    lambda_MLE        
    [1] 0.697252        
    lambda <- lambda_MLE
  3. Segunda iteração

    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

Agora itere até a convergência:

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

Nossa estimativa populacional é 1003774 e nossa taxa de poisson é estimada em 0,1994473 - essa é a proporção estimada da população amostrada. O principal problema que você terá nos problemas biológicos típicos com os quais está lidando é a suposição de que a taxa de poisson é uma constante.

Desculpe pelo post longo - este wiki não é realmente adequado para o código R.

Thylacoleo
fonte
3
Realce seu código e clique no botão que se parece com números binários ...
Shane
8

Isso soa como uma forma de 'marcar e recapturar', também conhecida como 'captura-recaptura', uma técnica bem conhecida em ecologia (e alguns outros campos, como epidemiologia). Não é minha área, mas o artigo da Wikipedia sobre marca e recaptura parece razoável, embora sua situação não seja a que o método Lincoln-Petersen explicado lá se aplica.

Eu acho que o shabbychef é o caminho certo para a sua situação, mas usar a distribuição Poisson para aproximar o binômio provavelmente tornaria as coisas um pouco mais simples e deve ser uma aproximação muito boa se o tamanho da população for muito grande, como nos seus exemplos. Eu acho que obter uma expressão explícita para a estimativa de probabilidade máxima do tamanho da população deve ser bem simples (veja, por exemplo, a Wikipedia novamente ), embora eu não tenha tempo para descobrir os detalhes agora.

uma parada
fonte
5

nkkP=1kmmn(nm)Pm(1-P)n-mnnkm(1-P)1

PmmPm/Pm+1(k-1)m+1n-mk

shabbychef
fonte