Оцените размер популяции по количеству повторных наблюдений

13

Скажем, у меня 50 миллионов уникальных вещей, и я беру 10 миллионов образцов (с заменой) ... Первый прикрепленный график показывает, сколько раз я выбираю одну и ту же "вещь", что относительно редко население больше, чем моя выборка.

Однако, если моя популяция составляет всего 10 миллионов штук, и я беру 10 миллионов выборок, как показывает второй график, я чаще буду повторять одну и ту же вещь несколько раз.

У меня вопрос - из моей таблицы частотных наблюдений (данные в гистограммах) можно ли получить оценку исходного размера популяции, когда она неизвестна? И было бы здорово, если бы вы могли предоставить указатель на то, как это сделать в R.

альтернативный текст

Аарон Стэтхэм
источник
Смотрите space.stackexchange.com/questions/41547/… для интересного приложения
kjetil b halvorsen

Ответы:

10

Как Гарван?

Проблема в том, что мы не знаем, сколько отсчетов наблюдается. Мы должны оценить это. Классической статистической процедурой для подобных ситуаций является алгоритм максимизации ожидания.

Простой пример:

Предположим, мы берем неизвестную популяцию (из 1000000) с постоянной Пуассона 0,2.

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

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

Но мы не наблюдаем нулевые счета. Вместо этого мы наблюдаем это:

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

table

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

Наблюдаются возможные частоты

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

Инициализируйте среднее значение распределения Пуассона - просто догадайтесь (мы знаем, что здесь 0,2).

lambda <- 1 
  1. Ожидание - распределение Пуассона

    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. максимизация

    lambda_MLE <- (1/sum(table))*(sum(table*k))        
    lambda_MLE        
    [1] 0.697252        
    lambda <- lambda_MLE
  3. Вторая итерация

    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

Теперь итерации до сходимости:

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

Наша оценка численности населения составляет 1003774, а наш коэффициент пуассона оценивается в 0,1994473 - это приблизительная доля выборки населения. Основная проблема, с которой вы столкнетесь в типичных биологических проблемах, с которыми вы сталкиваетесь, заключается в предположении, что скорость Пуассона является постоянной

Извините за скучный пост - эта вики не очень подходит для кода R.

Thylacoleo
источник
3
Выделите свой код и нажмите на кнопку, которая выглядит как двоичные числа ...
Шейн
8

Это звучит как форма «отлова и повторного захвата», известная как метод экологии (и некоторые другие области, такие как эпидемиология). Не моя сфера, а статья в Википедии о марке и возвращении выглядит разумной, хотя ваша ситуация не относится к той, к которой применяется метод Линкольна-Петерсена.

Я думаю, что shabbychef - это правильный путь для вашей ситуации, но использование распределения Пуассона для аппроксимации бинома, вероятно, сделает вещи немного проще и должно быть очень хорошим приближением, если численность населения очень велика, как в ваших примерах. Я думаю, что получение точного выражения для оценки максимального правдоподобия численности населения должно быть довольно простым (см., Например, снова Википедию ), хотя у меня нет времени, чтобы проработать детали прямо сейчас.

универсальный
источник
5

nkkP=1kmmn(nm)Pm(1P)nmnnkm(1P)1

PmmPm/Pm+1(k1)m+1nmk

shabbychef
источник