Saya memiliki beberapa frekuensi permintaan, dan saya perlu memperkirakan koefisien hukum Zipf. Ini adalah frekuensi teratas:
26486
12053
5052
3033
2536
2391
1444
1220
1152
1039
Saya memiliki beberapa frekuensi permintaan, dan saya perlu memperkirakan koefisien hukum Zipf. Ini adalah frekuensi teratas:
26486
12053
5052
3033
2536
2391
1444
1220
1152
1039
Jawaban:
Memperbarui Saya telah memperbarui kode dengan estimator kemungkinan maksimum sesuai saran @whuber. Meminimalkan jumlah kuadrat perbedaan antara probabilitas teoretis log dan frekuensi log meskipun memberikan jawaban akan menjadi prosedur statistik jika dapat ditunjukkan bahwa itu adalah semacam M-estimator. Sayangnya saya tidak bisa memikirkan yang bisa memberikan hasil yang sama.
Ini usahaku. Saya menghitung logaritma frekuensi dan mencoba menyesuaikannya dengan logaritma probabilitas teoretis yang diberikan oleh rumus ini . Hasil akhirnya tampak masuk akal. Ini kode saya di R.
Fit kuadratik terbaik adalah .s=1.47
Kemungkinan maksimum dalam R dapat dilakukan dengan
mle
fungsi (daristats4
paket), yang membantu menghitung kesalahan standar (jika fungsi kemungkinan maksimum negatif yang benar diberikan):Berikut adalah grafik kecocokan dalam skala log-log (lagi seperti yang disarankan @whuber):
Garis merah adalah jumlah kotak kuadrat, garis hijau adalah fit maksimum-likelihood.
sumber
Ada beberapa masalah sebelum kita dalam setiap masalah estimasi:
Perkirakan parameter.
Nilai kualitas estimasi itu.
Jelajahi data.
Evaluasi kecocokan.
Bagi mereka yang akan menggunakan metode statistik untuk memahami dan berkomunikasi, yang pertama tidak boleh dilakukan tanpa yang lain.
Dengan demikian probabilitas log untuk data adalah
Mengingat sifat hukum Zipf, cara yang tepat untuk membuat grafik kecocokan ini adalah pada plot log-log , di mana kecocokannya akan linear (menurut definisi):
Untuk mengevaluasi kebaikan kecocokan dan mengeksplorasi data, lihat residu (data / kecocokan, log-log sumbu lagi):
Karena residu tampak acak, dalam beberapa aplikasi kami mungkin puas untuk menerima Hukum Zipf (dan perkiraan parameter kami) sebagai deskripsi walaupun frekuensi kasar dapat diterima . Namun, analisis ini menunjukkan bahwa akan menjadi kesalahan untuk menganggap bahwa perkiraan ini memiliki nilai penjelas atau prediksi untuk set data yang diperiksa di sini.
sumber
Perkiraan Kemungkinan Maksimum hanya perkiraan titik parameters . Upaya ekstra diperlukan untuk menemukan juga interval kepercayaan dari estimasi tersebut. Masalahnya adalah bahwa interval ini tidak probabilistik. Orang tidak dapat mengatakan "nilai parameter s = ... adalah dengan probabilitas 95% dalam kisaran [...]".
Salah satu bahasa pemrograman probabilistik seperti PyMC3 membuat estimasi ini relatif mudah. Bahasa lain termasuk Stan yang memiliki fitur hebat dan komunitas yang mendukung.
Berikut ini adalah implementasi model Python saya yang dipasang pada data OPs (juga pada Github ):
Berikut perkiraan parameters dalam bentuk distribusi. Perhatikan betapa kompaknya perkiraan tersebut! Dengan probabilitas 95% nilai sebenarnya dari parameters berada dalam kisaran [1,439,1.461]; rerata adalah sekitar 1,45, yang sangat dekat dengan perkiraan MLE.
Untuk memberikan beberapa diagnosa pengambilan sampel dasar, kita dapat melihat bahwa pengambilan sampel "berbaur dengan baik" karena kita tidak melihat struktur apa pun dalam jejak:
Untuk menjalankan kode, kita perlu Python dengan paket Theano dan PyMC3 diinstal.
Terima kasih kepada @ w-huber atas jawaban dan komentarnya yang luar biasa!
sumber
Berikut ini adalah upaya saya untuk mencocokkan data, mengevaluasi dan mengeksplorasi hasil menggunakan VGAM:
Dalam kasus kami hipotesis nol Chi square adalah bahwa data didistribusikan sesuai dengan hukum zipf, maka nilai-p yang lebih besar mendukung klaim bahwa data didistribusikan sesuai dengan itu. Perhatikan bahwa bahkan nilai p yang sangat besar bukanlah bukti, hanya sebuah indikator.
sumber
Hanya untuk bersenang-senang, ini adalah contoh lain di mana UWSE dapat memberikan solusi formulir tertutup hanya menggunakan frekuensi paling atas - meskipun dengan biaya akurasi. Probabilitas menyalax = 1 unik di seluruh nilai parameter. Jikawx = 1^ menunjukkan frekuensi relatif yang sesuai,
Dalam hal ini, sejakwx = 1^= 0.4695599775 , kita mendapatkan:
Sekali lagi, UWSE hanya menyediakan estimasi yang konsisten - tidak ada interval kepercayaan, dan kita dapat melihat beberapa trade-off dalam akurasi. solusi mpiktas di atas juga merupakan aplikasi dari UWSE - meskipun pemrograman diperlukan. Untuk penjelasan lengkap tentang penaksir, lihat: https://paradsp.wordpress.com/ - semuanya ada di bagian bawah.
sumber
Solusi saya mencoba untuk melengkapi jawaban yang diberikan oleh mpiktas dan whuber melakukan implementasi dengan Python. Frekuensi dan rentang x kami adalah:
Karena fungsi kita tidak didefinisikan dalam semua rentang, kita perlu memeriksa bahwa kita menormalkan setiap kali kita menghitungnya. Dalam kasus diskrit, pendekatan sederhana adalah dengan membagi dengan jumlah semua y (x). Dengan cara ini kita dapat membandingkan berbagai parameter.
Hasilnya memberi kami kemiringan 1.450408 seperti pada jawaban sebelumnya.
sumber