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,