Estimando o tamanho da população a partir da frequência de duplicados e únicos amostrados

Existe um serviço da web em que posso solicitar informações sobre um item aleatório. Para cada solicitação, cada item tem a mesma chance de ser devolvido.

Posso continuar solicitando itens e registrar o número de duplicatas e único. Como posso usar esses dados para estimar o número total de itens?

probability population coupon-collector-problem hoju
fonte

O que você deseja estimar não é um tamanho de amostra, mas o tamanho de uma população (número total de itens exclusivos retornados pelo serviço da web).

GaBorgulya

Respostas:

Essa é essencialmente uma variante do problema do coletor de cupons.

Se houver itens no total e você tiver retirado um tamanho de amostra com substituição, a probabilidade de identificar itens únicos é $n$ $s$ $u$ onde

P r (U = u | n, s) = \frac{S_{2} (s, u) n!}{(n - u)! n^{s}}

$Pr(U=u|n,s) = \frac{S_2(s,u) n! }{ (n-u)! n^s }$

fornecenúmeros Stirling do segundo tipo

S_{2} (s, u)

$S_2(s,u)$

Agora tudo o que você precisa é de uma distribuição prévia para , aplicar teorema de Bayes, e obter uma distribuição posterior para . $Pr(N=n)$ $N$

Henry
fonte

Isso parece perder algumas informações, porque não leva em consideração as frequências com as quais os itens foram observados 2, 3, 4, ... vezes.

whuber

@ whuber: Pode parecer não usar as informações, mas se você investigar mais, deve descobrir que o número de itens exclusivos é uma estatística suficiente. Por exemplo, se você tirar uma amostra com a substituição de 4 itens de uma população de

n

$n$ , a probabilidade de obter 3 de um item e 1 de outro é

a de obter 2 cada um dos dois itens, não importa o queseja

; portanto, conhecer as frequências detalhadas não fornece informações mais úteis sobre a população do que simplesmente saber que havia dois itens únicos encontrados na amostra.

\frac{4}{3}

$\frac{4}{3}$

n

$n$

Henry

Ponto interessante sobre a suficiência do número de itens exclusivos. Portanto, as frequências podem servir como uma verificação das suposições (de independência e probabilidade igual), mas, caso contrário, são desnecessárias.

whuber

Já dei uma sugestão baseada nos números de Stirling do segundo tipo e nos métodos bayesianos.

Para aqueles que consideram os números de Stirling muito grandes ou os métodos bayesianos muito difíceis, um método mais difícil pode ser o de usar

E [U | n, s] = n (1 - {(1 - \frac{1}{n})}^{s})

$E[U|n,s] = n\left( 1- \left(1-\frac{1}{n}\right)^s\right)$

v a r [U | n, s] = n {(1 - \frac{1}{n})}^{s} + n^{2} (1 - \frac{1}{n}) {(1 - \frac{2}{n})}^{s} - n^{2} {(1 - \frac{1}{n})}^{2 s}

$var[U|n,s] = n\left(1-\frac{1}{n}\right)^s + n^2 \left(1-\frac{1}{n}\right)\left(1-\frac{2}{n}\right)^s - n^2\left(1-\frac{1}{n}\right)^{2s}$

e calcular novamente usando métodos numéricos.

Por exemplo, tomando o exemplo de GaBorgulya com e um observado , isso pode nos dar uma estimativa do $s=300$ $U = 265$ $\hat{n} \approx 1180$

Se essa fosse a população, isso nos daria uma variação de de cerca de 25, e dois desvios-padrão arbitrários de ambos os lados de 265 seriam de 255 e 275 (como eu disse, esse é um método grosseiro). 255 nos daria uma estimativa para cerca de 895, enquanto 275 daria cerca de 1692. O exemplo 1000 é confortavelmente dentro desse intervalo. $U$ $n$

Henry
fonte

s / n

$s/n$

n

$n$

n

$n$

s / n

$s/n$

U

$U$

1 - (1 - 1 / n)^{s} \geq (1 - f_{k} (s / n)) / f_{k} (s / n)

$1 - (1-1/n)^s \geq (1-f_k(s/n)) / f_k(s/n)$

f_{k} (x) = \sum_{i = 0}^{k} x^{i} / i!

$f_k(x) = \sum_{i=0}^k x^i/i!$

k

$k$

e^{x}

$e^x$

k = 1

$k=1$

\tilde{n} = \frac{s}{s - U} U

$\tilde{n} = \frac{s}{s-U} U$

s

$s$

\hat{n}

$\hat{n}$

Você pode usar o método de captura-recaptura , também implementado como o pacote Rcapture R .

Aqui está um exemplo, codificado em R. Vamos supor que o serviço da Web tenha N = 1000 itens. Faremos n = 300 solicitações. Gere uma amostra aleatória em que, numerando os elementos de 1 a k, onde k é quantos itens diferentes vimos.

N = 1000; population = 1:N # create a population of the integers from 1 to 1000
n = 300 # number of requests
set.seed(20110406)
observation = as.numeric(factor(sample(population, size=n,
  replace=TRUE))) # a random sample from the population, renumbered
table(observation) # a table useful to see, not discussed
k = length(unique(observation)) # number of unique items seen
(t = table(table(observation)))

O resultado da simulação é

  1   2   3 
234  27   4

assim, entre as 300 solicitações, havia 4 itens vistos 3 vezes, 27 itens vistos duas vezes e 234 itens vistos apenas uma vez.

Agora estime N desta amostra:

require(Rcapture)
X = data.frame(t)
X[,1]=as.numeric(X[,1])
desc=descriptive(X, dfreq=TRUE, dtype="nbcap", t=300)
desc # useful to see, not discussed
plot(desc) # useful to see, not discussed
cp=closedp.0(X, dfreq=TRUE, dtype="nbcap", t=300, trace=TRUE)
cp

O resultado:

Number of captured units: 265 

Abundance estimations and model fits:
                  abundance       stderr      deviance   df           AIC
M0**                  265.0          0.0  2.297787e+39  298  2.297787e+39
Mh Chao              1262.7        232.5  7.840000e-01    9  5.984840e+02
Mh Poisson2**         265.0          0.0  2.977883e+38  297  2.977883e+38
Mh Darroch**          553.9         37.1  7.299900e+01  297  9.469900e+01
Mh Gamma3.5**  5644623606.6  375581044.0  5.821861e+05  297  5.822078e+05

 ** : The M0 model did not converge
 ** : The Mh Poisson2 model did not converge
 ** : The Mh Darroch model did not converge
 ** : The Mh Gamma3.5 model did not converge
Note: 9 eta parameters has been set to zero in the Mh Chao model

$\hat{N}$

EDIT: Para verificar a confiabilidade do método acima, executei o código acima em 10000 amostras geradas. O modelo Mh Chao convergia todas as vezes. Aqui está o resumo:

> round(quantile(Nhat, c(0, 0.025, 0.25, 0.50, 0.75, 0.975, 1)), 1)
    0%   2.5%    25%    50%    75%  97.5%   100% 
 657.2  794.6  941.1 1034.0 1144.8 1445.2 2162.0 
> mean(Nhat)
[1] 1055.855
> sd(Nhat)
[1] 166.8352

GaBorgulya
fonte

Parece que é necessária alguma justificativa para o uso de modelos de captura e recuperação, porque esse não é um experimento padrão de captura e recuperação. (Possivelmente ele pode ser visto como 300 eventos de captura, mas a chamada para closedp não parecem indicar que.)

whuber

@whuber Sim, eu vi o exemplo como 300 eventos de captura. Como você quer dizer que "a chamada para closedp parece não indicar isso"? Aprecio as críticas (construtivas) e fico feliz em corrigir (ou excluir, se necessário) minha resposta, se ela estiver errada.

GaBorgulya

isso parece uma abordagem razoável. No entanto, não vou usar R, por isso preciso entender a matemática por trás disso. A página wiki cobre uma situação de 2 eventos - como eu a aplico neste caso?

Hoju

@ Ga vejo: você criou uma matriz de 300 x 300 para os dados! A ineficiência desse código me enganou: seria mais simples e mais direto usar `closedp.0 (Y, dfreq = TRUE, dtype =" nbcap ", t = 300) 'onde Y é a matriz de frequência {{1.234}, {2,27}, {3,4}} (que você calculou duas vezes e realmente exibiu!). Mais exatamente, as falhas de convergência são alarmantes, sugerindo que há problemas com o código ou os modelos subjacentes. (Uma busca exaustiva nos documentos por "M0" não mostra referências ou descrição para esse método ...)

whuber

@ whuber Eu simplifiquei o código seguindo a sua sugestão (dfreq = TRUE, dtype = "nbcap", t = 300). Obrigado novamente.

GaBorgulya