Perkirakan ukuran populasi yang dijadikan sampel dengan jumlah pengamatan berulang

13

Katakanlah saya memiliki populasi 50 juta hal unik, dan saya mengambil 10 juta sampel (dengan penggantian) ... Grafik pertama yang saya lampirkan menunjukkan berapa kali saya mencicipi "benda" yang sama, yang relatif jarang seperti populasi lebih besar dari sampel saya.

Namun jika populasi saya hanya 10 juta hal, dan saya mengambil 10 juta sampel, seperti grafik kedua menunjukkan saya akan lebih sering mencicipi hal yang sama berulang kali.

Pertanyaan saya adalah - dari tabel frekuensi pengamatan saya (data dalam diagram batang) apakah mungkin untuk mendapatkan estimasi ukuran populasi asli ketika tidak diketahui? Dan akan lebih bagus jika Anda bisa memberikan petunjuk bagaimana cara melakukannya di R.

teks alternatif

Aaron Statham
sumber

Jawaban:

10

Bagaimana kabar Garvan?

Masalahnya adalah kita tidak tahu berapa jumlah nol yang diamati. Kami harus memperkirakan ini. Prosedur statistik klasik untuk situasi seperti ini adalah algoritma Ekspektasi-Maksimalisasi.

Contoh sederhana:

Asumsikan kita menarik dari populasi yang tidak diketahui (1.000.000) dengan konstanta poisson 0,2.

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

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

Tapi kami tidak mengamati jumlah nol. Sebaliknya, kami mengamati ini:

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

table

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

Frekuensi yang mungkin diamati

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

Inisialisasi mean distribusi Poisson - cukup tebak (kami tahu 0,2 di sini).

lambda <- 1 
  1. Harapan - Distribusi 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. Maksimalisasi

    lambda_MLE <- (1/sum(table))*(sum(table*k))        
    lambda_MLE        
    [1] 0.697252        
    lambda <- lambda_MLE
  3. Iterasi kedua

    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

Sekarang beralih hingga konvergensi:

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

Perkiraan populasi kami adalah 1003774 dan laju poisson kami diperkirakan sebesar 0,1994473 - ini adalah proporsi estimasi populasi sampel. Masalah utama yang akan Anda miliki dalam masalah biologis tipikal yang Anda hadapi adalah asumsi bahwa tingkat poisson adalah konstan.

Maaf untuk posting yang bertele-tele - wiki ini tidak benar-benar cocok untuk kode R.

Thylacoleo
sumber
3
Sorot kode Anda dan klik tombol yang terlihat seperti angka biner ...
Shane
8

Ini terdengar seperti bentuk 'tanda dan tangkap kembali' alias 'tangkap-tangkap kembali', teknik yang terkenal di bidang ekologi (dan beberapa bidang lain seperti epidemiologi). Bukan wilayah saya, tetapi artikel Wikipedia tentang tanda dan penangkapan kembali terlihat masuk akal, meskipun situasi Anda bukanlah situasi yang dijelaskan oleh metode Lincoln – Petersen di sana.

Saya pikir shabbychef adalah salah satu jalur yang tepat untuk situasi Anda, tetapi menggunakan distribusi Poisson untuk memperkirakan binomial mungkin akan membuat segalanya lebih sederhana dan harus menjadi perkiraan yang sangat baik jika ukuran populasi sangat besar, seperti pada contoh Anda. Saya pikir mendapatkan ekspresi eksplisit untuk perkiraan kemungkinan maksimum ukuran populasi harus cukup mudah (lihat misalnya Wikipedia lagi ), meskipun saya tidak punya waktu untuk menyelesaikan rinciannya sekarang.

onestop
sumber
5

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

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

shabbychef
sumber