Bagaimana cara menghitung koefisien hukum Zipf dari satu set frekuensi teratas?

25

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
Diegolo
sumber
menurut halaman wikipedia , hukum Zipf memiliki dua parameter. Jumlah elemen dan s eksponen. Apa N dalam kasus Anda, 10? Dan frekuensi dapat dihitung dengan membagi nilai yang Anda berikan dengan jumlah semua nilai yang disediakan? NsN
mpiktas
biarkan sepuluh, dan frekuensi dapat dihitung dengan membagi nilai yang Anda berikan dengan jumlah semua nilai yang disediakan .. bagaimana saya bisa memperkirakan?
Diegolo

Jawaban:

22

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.

fr <- c(26486, 12053, 5052, 3033, 2536, 2391, 1444, 1220, 1152, 1039)

p <- fr/sum(fr)

lzipf <- function(s,N) -s*log(1:N)-log(sum(1/(1:N)^s))

opt.f <- function(s) sum((log(p)-lzipf(s,length(p)))^2)

opt <- optimize(opt.f,c(0.5,10))

> opt
$minimum
[1] 1.463946

$objective
[1] 0.1346248

Fit kuadratik terbaik adalah .s=1.47

Kemungkinan maksimum dalam R dapat dilakukan dengan mlefungsi (dari stats4paket), yang membantu menghitung kesalahan standar (jika fungsi kemungkinan maksimum negatif yang benar diberikan):

ll <- function(s) sum(fr*(s*log(1:10)+log(sum(1/(1:10)^s))))

fit <- mle(ll,start=list(s=1))

> summary(fit)
Maximum likelihood estimation

Call:
mle(minuslogl = ll, start = list(s = 1))

Coefficients:
  Estimate  Std. Error
s 1.451385 0.005715046

-2 log L: 188093.4 

Berikut adalah grafik kecocokan dalam skala log-log (lagi seperti yang disarankan @whuber):

s.sq <- opt$minimum
s.ll <- coef(fit)

plot(1:10,p,log="xy")
lines(1:10,exp(lzipf(s.sq,10)),col=2)
lines(1:10,exp(lzipf(s.ll,10)),col=3)

Garis merah adalah jumlah kotak kuadrat, garis hijau adalah fit maksimum-likelihood.

Log-log graph of fits

mpiktas
sumber
1
Ada juga paket R zipfR cran.r-project.org/web/packages/zipfR/index.html saya belum mencobanya.
onestop
@onestop, terima kasih atas tautannya. Akan lebih baik jika seseorang menjawab pertanyaan ini menggunakan paket ini. Solusi saya pasti tidak memiliki kedalaman, meskipun memberikan semacam jawaban.
mpiktas
(+1) Kamu sangat mengesankan. Begitu banyak kontribusi bagus di berbagai bidang statistik!
chl
@ chl, terima kasih! Saya tentu merasa bahwa saya bukan satu-satunya dengan karakteristik seperti itu di situs ini;)
mpiktas
25

Ada beberapa masalah sebelum kita dalam setiap masalah estimasi:

  1. Perkirakan parameter.

  2. Nilai kualitas estimasi itu.

  3. Jelajahi data.

  4. Evaluasi kecocokan.

Bagi mereka yang akan menggunakan metode statistik untuk memahami dan berkomunikasi, yang pertama tidak boleh dilakukan tanpa yang lain.

i=1,2,,nisss>0

Hs(n)=11s+12s++1ns.

i1n

log(Pr(i))=log(isHs(n))=slog(i)log(Hs(n)).

fi,i=1,2,,n

Pr(f1,f2,,fn)=Pr(1)f1Pr(2)f2Pr(n)fn.

Dengan demikian probabilitas log untuk data adalah

Λ(s)=si=1nfilog(i)(i=1nfi)log(Hs(n)).

s Kemungkinan .

s^=1.45041Λ(s^)=94046.7s^ls=1.463946Λ(s^ls)=94049.5

s[1.43922,1.46162] (jika saya melakukan perhitungan dengan benar :-).

Mengingat sifat hukum Zipf, cara yang tepat untuk membuat grafik kecocokan ini adalah pada plot log-log , di mana kecocokannya akan linear (menurut definisi):

enter image description here

Untuk mengevaluasi kebaikan kecocokan dan mengeksplorasi data, lihat residu (data / kecocokan, log-log sumbu lagi):

enter image description here

χ2=656.476 .


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.

whuber
sumber
1
@whuber, saya mungkin dengan rendah hati menyarankan sedikit hati-hati dengan formulasi yang diberikan di atas. Hukum Zipf biasanya dinyatakan sebagai hasil frekuensi relatif. Ini bukan (biasanya dianggap) distribusi dari mana sampel iid diambil. Kerangka kerja iid mungkin bukan ide terbaik untuk data ini. Mungkin saya akan memposting lebih lanjut tentang ini nanti.
kardinal
3
@ kardinal Saya menantikan apa yang Anda katakan. Jika Anda tidak punya waktu untuk tanggapan menyeluruh, bahkan sketsa apa yang menurut Anda mungkin merupakan "ide terbaik untuk data ini" akan sangat diterima. Saya bisa menebak di mana Anda akan pergi dengan ini: data telah diberi peringkat, sebuah proses yang menciptakan dependensi dan harus meminta saya untuk mempertahankan kemungkinan yang diperoleh tanpa mengenali efek potensial dari peringkat. Akan menyenangkan untuk melihat prosedur estimasi dengan justifikasi yang lebih baik. Saya berharap, bagaimanapun, bahwa analisis saya dapat diselamatkan oleh ukuran semata-mata dataset.
whuber
1
@ cardinal, jangan lakukan Fermat pada kami :) Jika Anda memiliki wawasan yang berbeda dari penjawab lainnya, jangan ragu untuk mengungkapkannya di jawaban yang terpisah, meskipun itu bukan merupakan jawaban yang sah saja. Dalam math.SE misalnya situasi seperti itu muncul cukup sering.
mpiktas
1
@ kardinal dengan mudah. Misalnya, Anda mengumpulkan frekuensi dan mengidentifikasi dan memberi peringkat sepuluh tertinggi. Anda berhipotesis tentang Hukum Zipf. Anda mengumpulkan satu set frekuensi baru dan melaporkannya berdasarkan peringkat sebelumnya . Itulah situasi awal, analisis saya sangat cocok, bergantung pada peringkat baru yang setuju dengan yang lama.
whuber
1
@whuber, terima kasih atas kesabaran Anda. Sekarang saya sepenuhnya mengerti alasan Anda. Di bawah model pengambilan sampel yang sekarang telah Anda sempurnakan sepenuhnya, saya setuju dengan analisis Anda. Mungkin pernyataan terakhir Anda masih agak licin. Jika penyortiran tidak menyebabkan ketergantungan yang kuat daripada metode Anda akan menjadi konservatif. Jika ketergantungan yang diinduksi cukup kuat, itu mungkin menjadi antikonservatif. Terima kasih atas kesabaran Anda di hadapan ilmu cabul saya.
kardinal
2

Perkiraan Kemungkinan Maksimum hanya perkiraan titik parameter s. 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 ):

import theano.tensor as tt
import numpy as np
import pymc3 as pm
import matplotlib.pyplot as plt

data = np.array( [26486, 12053, 5052, 3033, 2536, 2391, 1444, 1220, 1152, 1039] )

N = len( data )

print( "Number of data points: %d" % N )

def build_model():
    with pm.Model() as model:
        # unsure about the prior...
        #s = pm.Normal( 's', mu=0.0, sd=100 )
        #s = pm.HalfNormal( 's', sd=10 )
        s = pm.Gamma('s', alpha=1, beta=10)

        def logp( f ):
            r = tt.arange( 1, N+1 )
            return -s * tt.sum( f * tt.log(r) ) - tt.sum( f ) * tt.log( tt.sum(tt.power(1.0/r,s)) )

        pm.DensityDist( 'obs', logp=logp, observed={'f': data} )

    return model


def run( n_samples=10000 ):
    model = build_model()
    with model:
        start = pm.find_MAP()
        step = pm.NUTS( scaling=start )
        trace = pm.sample( n_samples, step=step, start=start )

    pm.summary( trace )
    pm.traceplot( trace )
    pm.plot_posterior( trace, kde_plot=True )
    plt.show()

if __name__ == '__main__':
    run()

Berikut perkiraan parameter sdalam bentuk distribusi. Perhatikan betapa kompaknya perkiraan tersebut! Dengan probabilitas 95% nilai sebenarnya dari parametersberada dalam kisaran [1,439,1.461]; rerata adalah sekitar 1,45, yang sangat dekat dengan perkiraan MLE.

enter image description here

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:

enter image description here

Untuk menjalankan kode, kita perlu Python dengan paket Theano dan PyMC3 diinstal.

Terima kasih kepada @ w-huber atas jawaban dan komentarnya yang luar biasa!

Vladislavs Dovgalecs
sumber
1

Berikut ini adalah upaya saya untuk mencocokkan data, mengevaluasi dan mengeksplorasi hasil menggunakan VGAM:

require("VGAM")

freq <- dzipf(1:100, N = 100, s = 1)*1000 #randomizing values
freq <- freq  + abs(rnorm(n=1,m=0, sd=100)) #adding noize

zdata <- data.frame(y = rank(-freq, ties.method = "first") , ofreq = freq)
fit = vglm(y ~ 1, zipf, zdata, trace = TRUE,weight = ofreq,crit = "coef")
summary(fit)

s <- (shat <- Coef(fit)) # the coefficient we've found
probs <- dzipf(zdata$y, N = length(freq), s = s) # expected values
chisq.test(zdata$ofreq, p = probs) 
plot(zdata$y,(zdata$ofreq),log="xy") #log log graph
lines(zdata$y, (probs)*sum(zdata$ofreq),  col="red") # red line, num of predicted frequency

enter image description here

    Chi-squared test for given probabilities

data:  zdata$ofreq
X-squared = 99.756, df = 99, p-value = 0.4598

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.

Guy s
sumber
0

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=1unik di seluruh nilai parameter. Jikawx=1^ menunjukkan frekuensi relatif yang sesuai,

sUWSE^=H10-1(1wx=1^)

Dalam hal ini, sejak wx=1^=0,4695599775, kita mendapatkan:

sUWSE^=1.4

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.

CYP450
sumber
Bagaimana hubungan UWSE dengan hukum Zipf?
Michael R. Chernick
UWSE (Unique Weight Space Estimation) menggunakan fakta bahwa probabilitas / frekuensi teratas adalah unik di berbagai nilai parameter s, untuk N yang diberikan, untuk menemukan s. Sehubungan dengan hukum Zipf, ini memberitahu kita bahwa dengan memberikan sejumlah item ke peringkat, N, dan frekuensi paling atas, hanya ada satu cara menetapkan frekuensi ke item yang tersisa (2, ..., N) sedemikian rupa sehingga kita dapat katakan "item ke-n adalah 1 / n ^ s kali lebih besar dari item yang paling sering, untuk beberapa s". Dengan kata lain, mengingat info ini, hanya ada satu cara bagi hukum Zipf untuk memegang - tentu saja, dengan asumsi bahwa hukum Zipf memang berlaku.
CYP450
0

Solusi saya mencoba untuk melengkapi jawaban yang diberikan oleh mpiktas dan whuber melakukan implementasi dengan Python. Frekuensi dan rentang x kami adalah:

freqs = np.asarray([26486, 12053, 5052, 3033, 2536, 2391, 1444, 1220, 1152, 1039])
x = np.asarray([1, 2, 3, 4, 5 ,6 ,7 ,8 ,9, 10])

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.

f,ax = plt.subplots()
ax.plot(x, f1, 'o')
ax.set_xscale("log")
ax.set_yscale("log")

def loglik(b):  
    # Power law function
    Probabilities = x**(-b)

    # Normalized
    Probabilities = Probabilities/Probabilities.sum()

    # Log Likelihoood
    Lvector = np.log(Probabilities)

    # Multiply the vector by frequencies
    Lvector = np.log(Probabilities) * freqs

    # LL is the sum
    L = Lvector.sum()

    # We want to maximize LogLikelihood or minimize (-1)*LogLikelihood
    return(-L)

s_best = minimize(loglik, [2])
print(s_best)
ax.plot(x, freqs[0]*x**-s_best.x)

enter image description here

Hasilnya memberi kami kemiringan 1.450408 seperti pada jawaban sebelumnya.

ivangtorre
sumber