Bagaimana cara memperkirakan parameter untuk distribusi terpotong Zipf dari sampel data?

Saya punya masalah dengan parameter estimasi untuk Zipf. Situasi saya adalah sebagai berikut:

Saya memiliki kumpulan sampel (diukur dari percobaan yang menghasilkan panggilan yang harus mengikuti distribusi Zipf). Saya harus menunjukkan bahwa generator ini benar-benar menghasilkan panggilan dengan distribusi zipf. Saya sudah membaca T&J ini. Bagaimana cara menghitung koefisien hukum Zipf dari satu set frekuensi teratas? tapi saya mencapai hasil yang buruk karena saya menggunakan distribusi terpotong. Misalnya jika saya menetapkan nilai "s" menjadi "0,9" untuk proses pembuatan, jika saya mencoba memperkirakan nilai "s" seperti yang ditulis dalam Q&A yang dilaporkan, saya memperoleh "s" sama dengan 0,2 ca. Saya pikir ini adalah karena saya menggunakan distribusi TRUNCATED (saya harus membatasi zipf dengan titik pemotongan, itu benar-terpotong).

Bagaimana saya bisa memperkirakan parameter dengan distribusi zipf terpotong?

distributions estimation pareto-distribution zipf Maurizio
sumber

untuk menjadi jelas, apa yang Anda benar memotong kanan? Distribusi nilai atau plot Zipf itu sendiri? Apakah Anda tahu titik pemotongan? Apakah pemotongan merupakan artefak data atau artefak pemrosesan data (misalnya, beberapa keputusan yang Anda atau eksperimen lakukan)? Setiap detail tambahan akan sangat membantu.

kardinal

@kardinal. (bagian 1/2) Terima kasih kardinal. Saya akan memberikan rincian lebih lanjut: Saya memiliki generator VoIP yang menghasilkan panggilan mengikuti Zipf (dan distribusi lainnya) untuk volume per pemanggil. Saya harus memverifikasi bahwa generator ini benar-benar mengikuti distribusi ini. Untuk Distribusi Zipf saya perlu mendefinisikan titik pemotongan (oleh karena itu diketahui dan mengacu pada distribusi nilai) yang merupakan jumlah maksimum panggilan yang dihasilkan oleh pengguna dan parameter skala. Khususnya dalam kasus saya, nilai ini sama dengan 500, yang menunjukkan bahwa satu pengguna dapat menghasilkan 500 panggilan maksimum.

Maurizio

(bagian 2/2) Parameter lain yang akan ditetapkan adalah parameter skala untuk Zipf yang menentukan penyebaran distribusi (nilai ini dalam kasus saya adalah 0,9). Saya memiliki semua parameter (ukuran sampel, frekuensi per pengguna, dll) tetapi saya harus memverifikasi bahwa dataset saya mengikuti distribusi zipf.

Maurizio

jadi Anda tampaknya menormalkan kembali distribusi dengan , karena untuk, apa yang saya anggap sebagai "terpotong Zipf", parameter penskalaan 0,9 tidak mungkin dilakukan . Jika Anda dapat menghasilkan banyak data ini dan Anda "hanya" memiliki 500 hasil yang mungkin, mengapa tidak menggunakan uji chi-square good-of-fit? Karena distribusi Anda memiliki ekor yang panjang, Anda mungkin perlu ukuran sampel yang cukup besar. Tapi, itu akan menjadi satu arah. Metode cepat-dan-kotor lainnya adalah untuk memeriksa bahwa Anda mendapatkan distribusi empiris yang tepat untuk nilai - nilai kecil dari jumlah panggilan.

\sum_{i = 1}^{500} i^{- 0.9}

$\sum_{i=1}^{500} i^{-0.9}$

kardinal

Jawaban:

Pembaruan : 7 Apr 2011 Jawaban ini semakin panjang dan mencakup banyak aspek masalah yang dihadapi. Namun, saya telah menolak, sejauh ini, memecahnya menjadi jawaban yang terpisah.

Saya telah menambahkan di bagian paling bawah diskusi tentang kinerja Pearson's untuk contoh ini. $\chi^2$

Bruce M. Hill menulis, mungkin, makalah "seminal" tentang estimasi dalam konteks seperti Zipf. Dia menulis beberapa makalah pada pertengahan 1970 tentang topik itu. Namun, "penaksir Hill" (seperti yang sekarang disebut) pada dasarnya bergantung pada statistik urutan maksimal sampel dan, tergantung pada jenis pemotongan yang ada, yang dapat membuat Anda dalam beberapa masalah.

Makalah utama adalah:

BM Hill, Sebuah pendekatan umum sederhana untuk menarik kesimpulan tentang ekor suatu distribusi , Ann. Stat. , 1975.

Jika data Anda benar-benar awalnya Zipf dan kemudian dipotong, maka korespondensi yang bagus antara distribusi derajat dan plot Zipf dapat dimanfaatkan untuk keuntungan Anda.

Secara khusus, distribusi derajat hanyalah distribusi empiris dari berapa kali setiap respons bilangan terlihat,

d_{i} = \frac{# {j : X_{j} = i}}{n} .

$d_i = \frac{\#\{j: X_j = i\}}{n} .$

Jika kita plot ini terhadap pada plot log-log, kita akan mendapatkan tren linier dengan kemiringan yang sesuai dengan koefisien penskalaan. $i$

Di sisi lain, jika kita memplot plot Zipf , di mana kita mengurutkan sampel dari yang terbesar ke yang terkecil dan kemudian memplot nilai-nilai terhadap peringkat mereka, kita mendapatkan tren linier yang berbeda dengan kemiringan yang berbeda . Namun lereng terkait.

Jika adalah koefisien hukum skala untuk distribusi Zipf, maka kemiringan dalam plot pertama adalah dan kemiringan dalam plot kedua adalah . Di bawah ini adalah contoh plot untuk dan . Panel kiri adalah distribusi derajat dan kemiringan garis merah adalah . Sisi kanan adalah plot Zipf, dengan garis merah yang ditumpangkan memiliki kemiringan . $\alpha$ $-\alpha$ $-1/(\alpha-1)$ $\alpha = 2$ $n = 10^6$ $-2$ $-1/(2-1) = -1$

Distribusi derajat (kiri) dan plot Zipf (kanan) untuk sampel iid dari distribusi Zipf.

Jadi, jika data Anda telah terpotong sehingga Anda tidak melihat nilai yang lebih besar dari beberapa ambang , tetapi data tersebut didistribusikan secara Zipf dan cukup besar, maka Anda dapat memperkirakan dari distribusi derajat . Pendekatan yang sangat sederhana adalah mencocokkan baris ke plot log-log dan menggunakan koefisien yang sesuai. $\tau$ $\tau$ $\alpha$

Jika data Anda terpotong sehingga Anda tidak melihat nilai - nilai kecil (misalnya, cara banyak penyaringan dilakukan untuk set data web yang besar), maka Anda dapat menggunakan plot Zipf untuk memperkirakan kemiringan pada skala log-log dan kemudian " mundur "eksponen penskalaan. Katakanlah perkiraan kemiringan Anda dari plot Zipf adalah . Kemudian, satu perkiraan sederhana dari koefisien scaling-law adalah $\hat{\beta}$

\hat{α} = 1 - \frac{1}{\hat{β}} .

$\hat{\alpha} = 1 - \frac{1}{\hat{\beta}} .$

@csgillespie memberikan satu makalah baru yang ditulis bersama oleh Mark Newman di Michigan mengenai topik ini. Dia sepertinya menerbitkan banyak artikel serupa tentang ini. Di bawah ini adalah satu lagi bersama dengan beberapa referensi lain yang mungkin menarik. Newman terkadang tidak melakukan hal yang paling masuk akal secara statistik, jadi berhati-hatilah.

MEJ Newman, hukum Power, distribusi Pareto dan hukum Zipf , Fisika Kontemporer 46, 2005, hlm. 323-351.

M. Mitzenmacher, Sejarah Singkat Model Generatif untuk Hukum Daya dan Distribusi Lognormal , Matematika Internet. , vol. 1, tidak. 2, 2003, hlm. 226-251.

K. Knight, Sebuah modifikasi sederhana dari estimator Hill dengan aplikasi untuk ketahanan dan pengurangan bias , 2010.

Adendum :

$R$ $10^5$

> x <- (1:500)^(-0.9)
> p <- x / sum(x)
> y <- sample(length(p), size=100000, repl=TRUE, prob=p)
> tab <- table(y)
> plot( 1:500, tab/sum(tab), log="xy", pch=20, 
        main="'Truncated' Zipf simulation (truncated at i=500)",
        xlab="Response", ylab="Probability" )
> lines(p, col="red", lwd=2)

Plot yang dihasilkan adalah

Zipf "Terpotong" (terpotong pada i = 500)

$i \leq 30$

Namun, dari sudut pandang praktis, plot semacam itu harus relatif menarik.

$\alpha = 2$ $n = 300\,000$ $x_{\mathrm{max}} = 500$

$\chi^2$

X^{2} = \sum_{i = 1}^{500} \frac{(O_{i} - E_{i})^{2}}{E_{i}}

$X^2 = \sum_{i=1}^{500} \frac{(O_i - E_i)^2}{E_i}$

O_{i}

$O_i$

i

$i$

E_{i} = n p_{i} = n i^{- α} / \sum_{j = 1}^{500} j^{- α}

$E_i = n p_i = n i^{-\alpha} / \sum_{j=1}^{500} j^{-\alpha}$

Kami juga akan menghitung statistik kedua yang dibentuk dengan terlebih dahulu menampar hitungan dalam nampan berukuran 40, seperti yang ditunjukkan dalam lembar kerja Maurizio (nampan terakhir hanya berisi jumlah dari dua puluh nilai hasil terpisah.

$n$ $p$

$p$

masukkan deskripsi gambar di sini

$R$

# Chi-square testing of the truncated Zipf.

a <- 2
n <- 300000
xmax <- 500

nreps <- 5000

zipf.chisq.test <- function(n, a=0.9, xmax=500, bin.size = 40)
{
  # Make the probability vector
  x <- (1:xmax)^(-a)
  p <- x / sum(x)

  # Do the sampling
  y <- sample(length(p), size=n, repl=TRUE, prob=p)

  # Use tabulate, NOT table!
  tab <- tabulate(y,xmax)

  # unbinned chi-square stat and p-value
  discrepancy <- (tab-n*p)^2/(n*p)
  chi.stat <- sum(discrepancy)
  p.val    <- pchisq(chi.stat, df=xmax-1, lower.tail = FALSE)

  # binned chi-square stat and p-value
  bins <- seq(bin.size,xmax,by=bin.size)
  if( bins[length(bins)] != xmax )
    bins <- c(bins, xmax)

  tab.bin  <- cumsum(tab)[bins]
  tab.bin <- c(tab.bin[1], diff(tab.bin))

  prob.bin <- cumsum(p)[bins] 
  prob.bin <- c(prob.bin[1], diff(prob.bin))

  disc.bin <- (tab.bin - n*prob.bin)^2/(n * prob.bin)
  chi.stat.bin <- sum(disc.bin)
  p.val.bin <- pchisq(chi.stat.bin, df=length(tab.bin)-1, lower.tail = FALSE)

  # Return the binned and unbineed p-values
  c(p.val, p.val.bin, chi.stat, chi.stat.bin)
}

set.seed( .Random.seed[2] )

all <- replicate(nreps, zipf.chisq.test(n, a, xmax))

par(mfrow=c(2,1))
hist( all[1,], breaks=20, col="darkgrey", border="white",
      main="Histogram of unbinned chi-square p-values", xlab="p-value")
hist( all[2,], breaks=20, col="darkgrey", border="white",
      main="Histogram of binned chi-square p-values", xlab="p-value" )

type.one.error <- rowMeans( all[1:2,] < 0.05 )

kardinal
sumber

+1, jawaban yang bagus seperti biasa. Anda harus mencalonkan diri sebagai moderator, masih ada 1 jam lagi :)

mpiktas

@mpikta, terima kasih atas pujian dan dukungannya. Saya tidak yakin bisa membenarkan mencalonkan diri dengan kandidat yang sudah sangat kuat yang, secara seragam, berpartisipasi lebih luas dan lebih lama daripada saya.

kardinal

@ cardinal, berikut adalah beberapa tautan ke alternatif untuk penaksir Hill: artikel asli oleh Paulauskas dan tindak lanjut oleh Vaiciulis dan Gadeikis dan Paulauskas . Pengukur ini diduga memiliki sifat yang lebih baik daripada Hill asli.

mpiktas

@mpikta, terima kasih atas tautannya. Ada beberapa versi "penaksir Hill" yang baru dan lebih baik. Kelemahan utama dari pendekatan asli adalah bahwa ia membutuhkan pilihan "cutoff" di mana berhenti rata-rata. Saya pikir sebagian besar yang telah dilakukan oleh "eyeballing" itu yang membuka satu hingga tuduhan subjektivitas. Salah satu buku Resnick tentang distribusi berekor panjang membahas hal ini secara terperinci, jika saya ingat. Saya pikir ini yang terbaru.

kardinal

@ kardinal, terima kasih banyak, Anda sangat baik dan sangat detail! Contoh Anda dalam R sangat berguna bagi saya, tetapi bagaimana saya bisa melakukan uji chi-square formal dalam kasus ini? (Saya menggunakan uji chi-square dengan distribusi lain seperti seragam, eksponensial, normal, tetapi saya memiliki banyak keraguan tentang zipf..Maaf tapi ini adalah pendekatan pertama saya untuk topik ini). Pertanyaan untuk modetator: apakah saya harus menulis T&J lain seperti "bagaimana melakukan uji chi-square untuk distribusi zipf terpotong?" atau lanjutkan dalam Tanya Jawab ini, mungkin memperbarui tag dan judul?

Maurizio

Kertas

Clauset, A et al , Distribusi Power-law dalam Data Empiris . 2009

berisi deskripsi yang sangat baik tentang bagaimana cara menyesuaikan model hukum kekuasaan Halaman web terkait memiliki sampel kode. Sayangnya, itu tidak memberikan kode untuk distribusi terpotong, tetapi mungkin memberi Anda pointer.

Sebagai tambahan, makalah ini membahas fakta bahwa banyak "dataset kekuasaan-hukum" dapat dimodelkan dengan baik (dan dalam beberapa kasus lebih baik) dengan distribusi Log normal atau eksponensial!

csgillespie
sumber

Sayangnya makalah ini tidak mengatakan apa-apa tentang distribusi terpotong .. Saya menemukan beberapa paket dalam R yang berurusan dengan parameter estimasi Zipf dengan cara yang sederhana (zipfR, VGAM) tetapi distribusi terpotong memerlukan "perlakuan khusus". Dengan kalimat terakhir Anda, apakah yang Anda maksudkan adalah mungkin untuk membuat model dataset hukum-kekuatan dengan distribusi eksponensial misalnya dan kemudian menerapkan beberapa proses parameter estimasi untuk distribusi eksponensial "terpotong"? Saya sangat pemula dalam topik ini!

Maurizio

Dalam makalah, penulis menganalisis ulang set data yang berbeda di mana kekuatan-hukum telah dipasang. Para penulis menunjukkan bahwa dalam sejumlah kasus model hukum-kekuasaan tidak terlalu bagus dan distribusi alternatif akan lebih baik.

csgillespie

Mengikuti jawaban terperinci dari kardinal pengguna, saya melakukan uji chi-square pada distribusi zipf saya yang mungkin terpotong. Hasil uji chi-square dilaporkan dalam tabel berikut:

masukkan deskripsi gambar di sini

Di mana StartInterval dan EndInterval mewakili misalnya rentang panggilan dan Yang Teramati adalah jumlah penelepon yang menghasilkan dari 0 hingga 19 panggilan, dan seterusnya .. Uji chi-square baik sampai kolom terakhir tercapai, mereka meningkatkan final perhitungan, jika tidak sampai titik itu nilai chi-square "parsial" dapat diterima!

Dengan tes lain hasilnya sama, kolom terakhir (atau 2 kolom terakhir) selalu meningkatkan nilai akhir dan saya tidak tahu mengapa dan saya tidak tahu jika (dan bagaimana) menggunakan tes validasi lain.

PS: untuk kelengkapan, untuk menghitung nilai yang diharapkan ( Diharapkan ) saya mengikuti saran kardinal dengan cara ini:

masukkan deskripsi gambar di sini

di mana x_i 's digunakan untuk menghitung: x <- (1:n)^-S, yang P_i ' s untuk menghitung p <- x / sum(x)dan akhirnya E_i (diharapkan nr pengguna untuk setiap nr panggilan) diperoleh denganP_i * Total_Caller_Observed

dan dengan Derajat Kebebasan = 13 kebaikan Chi-Square selalu menolak Hipotesis bahwa set sampel mengikuti Distribusi Zipf karena Statistik Uji (64,14 dalam kasus ini) lebih besar daripada yang dilaporkan dalam tabel chi-square, "demerit" untuk kolom terakhir. Hasil grafis dilaporkan di sini: masukkan deskripsi gambar di sini

meskipun titik pemotongan diatur ke 500 nilai maksimum yang diperoleh adalah 294. Saya pikir bahwa "dispersi" akhir adalah penyebab kegagalan uji chi-square.

MEMPERBARUI!!

Saya mencoba untuk melakukan uji chi-square pada sampel data zipf yang mungkin dihasilkan dengan kode R yang dilaporkan dalam jawaban di atas.

> x <- (1:500)^(-2)
> p <- x / sum(x)
> y <- sample(length(p), size=300000, repl=TRUE, prob=p)
> tab <- table(y)
> length(tab)
[1] 438
> plot( 1:438, tab/sum(tab), log="xy", pch=20, 
        main="'Truncated' Zipf simulation (truncated at i=500)",
        xlab="Response", ylab="Probability" )
> lines(p, col="red", lwd=2)

Plot terkait adalah sebagai berikut: masukkan deskripsi gambar di sini

Hasil uji chi-square dilaporkan dalam gambar berikut: masukkan deskripsi gambar di sini

dan statistik uji chi-square (44,57) terlalu tinggi untuk validasi dengan Gelar Kebebasan yang dipilih. Juga dalam hal ini "dispersi" data terakhir adalah penyebab dari nilai chi-square yang tinggi. Tetapi ada prosedur untuk memvalidasi distribusi zipf ini (terlepas dari generator "salah" saya, saya ingin fokus pada sampel data R) ???

Maurizio
sumber

@Maurizio, entah kenapa, saya melewatkan posting ini sampai sekarang. Apakah di sana Anda dapat mengeditnya dan menambahkan plot yang mirip dengan yang terakhir di posting saya, tetapi menggunakan data yang Anda amati? Itu mungkin membantu mendiagnosis masalah. Saya pikir saya melihat pertanyaan Anda yang lain di mana Anda mengalami masalah dalam menghasilkan distribusi yang seragam, jadi mungkin itu termasuk dalam analisis ini juga. (?) Salam.

kardinal

@ cardinal, saya memperbarui hasilnya! Bagaimana menurut anda? Pertanyaan tentang distribusi seragam adalah hal lain yang harus saya tentukan dengan cara yang lebih baik dan saya akan melakukannya hari ini atau besok;)

Maurizio

S = 0.9

$S = 0.9$

p = P (X_{i} = 500) \approx 4.05 \times 10^{- 4}

$p = \mathbb{P}(X_i = 500) \approx 4.05 \times 10^{-4}$

n = 8454

$n = 8454$

8454 \cdot 4.05 \cdot 10^{- 4} \approx 3.43

$8454 \cdot 4.05 \cdot 10^{-4} \approx 3.43$

1 - (1 - 0.000405)^{8454} \approx 0.9675

$1 - (1 - 0.000405)^{8454} \approx 0.9675$ . Perhatikan seberapa dekat yang cocok dengan simulasi di atas.

kardinal

@ cardinal, saya juga berpikir ada sesuatu yang "salah" dalam prosedur pembuatan (tujuan saya adalah untuk memvalidasi bahwa generator ini benar-benar mengikuti distribusi Zipf). Saya harus berbicara dengan para perancang proyek pada hari-hari ini.

Maurizio