Szacowanie wielkości populacji na podstawie częstotliwości duplikatów i niepowtarzalnych próbek

Istnieje serwis internetowy, w którym mogę poprosić o informacje na temat losowego przedmiotu. Przy każdym zamówieniu każdy przedmiot ma równą szansę na zwrot.

Potrafię nadal zamawiać przedmioty i rejestrować liczbę duplikatów i unikatowe. Jak mogę wykorzystać te dane do oszacowania całkowitej liczby produktów?

probability population coupon-collector-problem hoju
źródło

To, co chcesz oszacować, to nie wielkość próbki, ale wielkość populacji (całkowita liczba unikalnych elementów zwróconych przez sieć serc).

GaBorgulya

Odpowiedzi:

Jest to zasadniczo wariant problemu kolektora kuponów.

Jeżeli istnieje pozycji łącznie i miały wielkość próbki z wymianą to prawdopodobieństwo po zidentyfikowaniu unikalne elementów jest $n$ $s$ $u$ gdziedajeliczby Stirlinga drugiego rodzaju

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

S_{2} (s, u)

$S_2(s,u)$

Teraz wszystko, czego potrzebujemy, to przed dystrybucji , stosuje Twierdzenie Bayesa i dostać tylną dystrybucję dla . $Pr(N=n)$ $N$

Henz
źródło

Wydaje się, że traci to trochę informacji, ponieważ nie uwzględnia częstotliwości, z którymi obserwowano pozycje 2, 3, 4, ... razy.

whuber

@ whuber: Może się wydawać, że nie wykorzystujesz informacji, ale jeśli będziesz dalej badać, powinieneś stwierdzić, że liczba unikalnych przedmiotów jest wystarczającą statystyką. Na przykład, jeśli weźmiesz próbkę z wymianą 4 elementów z populacji

, prawdopodobieństwo uzyskania 3 jednego elementu i 1 innego wynosi

n

$n$

że otrzymanie 2 każdego z dwóch elementów, bez względu na to, czymjest

, więc znajomość szczegółowych częstotliwości nie daje więcej użytecznych informacji o populacji niż po prostu wiedza, że w próbie znaleziono dwa unikalne elementy.

\frac{4}{3}

$\frac{4}{3}$

n

$n$

Henry

Interesujący punkt dotyczący wystarczalności liczby unikalnych przedmiotów. Częstotliwości mogą więc służyć jako kontrola założeń (niezależności i równego prawdopodobieństwa), ale poza tym są niepotrzebne.

whuber

Dałem już sugestię opartą na liczbach Stirlinga drugiego rodzaju i metodach bayesowskich.

Dla tych, którzy uważają, że liczby Stirlinga są zbyt duże lub metody Bayesa zbyt trudne, może być zastosowanie bardziej surowej metody

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

i ponownie obliczyć za pomocą metod numerycznych.

$s=300$ $U = 265$ $\hat{n} \approx 1180$

$U$ $n$

Henz
źródło

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

Można użyć metody wychwytywania odbić , również wdrożony jako opakowania Rcapture R .

Oto przykład zakodowany w R. Załóżmy, że usługa sieciowa ma N = 1000 pozycji. Zrobimy n = 300 wniosków. Wygeneruj losową próbkę, gdzie, numerując elementy od 1 do k, gdzie k to ile różnych przedmiotów widzieliśmy.

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

Wynikiem symulacji jest

  1   2   3 
234  27   4

tak więc wśród 300 próśb były 4 pozycje widoczne 3 razy, 27 pozycji widziałem dwa razy, a 234 pozycji pokazano tylko raz.

Teraz oszacuj N z tej próbki:

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

Wynik:

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

EDYCJA: Aby sprawdzić wiarygodność powyższej metody, uruchomiłem powyższy kod na 10000 wygenerowanych próbkach. Model MH Chao zbierał się za każdym razem. Oto podsumowanie:

> 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
źródło

Wydaje się, że konieczne jest uzasadnienie zastosowania modeli przechwytywania i przechwytywania, ponieważ nie jest to standardowy eksperyment przechwytywania i przechwytywania. (Być może może to być postrzegane jako 300 zdarzeń przechwytywania, ale wezwanie do closedp nie wydaje się tego wskazywać).

whuber

@ whuber Tak, widziałem przykład 300 zdarzeń przechwytywania. Co masz na myśli mówiąc, że „wezwanie do closedp nie wydaje się tego wskazywać”? Doceniam (konstruktywną) krytykę i cieszę się, że mogę poprawić (lub usunąć, jeśli to konieczne) moją odpowiedź, jeśli okaże się ona błędna.

GaBorgulya,

wydaje się to rozsądnym podejściem. Jednak nie będę używać R, więc muszę zrozumieć matematykę. Strona wiki dotyczy sytuacji z 2 zdarzeniami - jak zastosować ją w tym przypadku?

hoju

@ Ga Widzę: Utworzyłeś macierz 300 x 300 dla danych! Nieefektywność tego kodu oszukała mnie: prostsze i bardziej bezpośrednie byłoby użycie `closedp.0 (Y, dfreq = TRUE, dtype =" nbcap ", t = 300) ', gdzie Y jest macierzą częstotliwości {{1,234}, {2,27}, {3,4}} (które obliczyłeś dwukrotnie i faktycznie wyświetliłeś!). Co więcej, niepowodzenia konwergencji są alarmujące, co sugeruje, że występują problemy z podstawowym kodem lub modelami. (Wyczerpujące przeszukanie dokumentacji dla „M0” nie

ujawniło

@ whuber Uprościłem kod zgodnie z twoją sugestią (dfreq = TRUE, dtype = "nbcap", t = 300). Dzięki jeszcze raz.

GaBorgulya,