Jak obliczyć współczynnik prawa Zipfa z zestawu najwyższych częstotliwości?

25

Mam kilka częstotliwości zapytań i muszę oszacować współczynnik prawa Zipfa. Oto najwyższe częstotliwości:

26486
12053
5052
3033
2536
2391
1444
1220
1152
1039
Diegolo
źródło
według strony wikipedii prawo Zipfa ma dwa parametry. Liczba elementów i wykładnik. Co to jest w twoim przypadku, 10? Czy częstotliwości można obliczyć, dzieląc dostarczone wartości przez sumę wszystkich podanych wartości? s NN.sN.
mpiktas
niech będzie dziesięć, a częstotliwości można obliczyć, dzieląc dostarczone wartości przez sumę wszystkich dostarczonych wartości .. jak mogę oszacować?
Diegolo

Odpowiedzi:

22

Aktualizacja Zaktualizowałem kod z estymatorem maksymalnego prawdopodobieństwa zgodnie z sugestią @whuber. Minimalizacja sumy kwadratów różnic między logicznymi prawdopodobieństwami teoretycznymi a częstotliwościami logów daje odpowiedź, byłaby to procedura statystyczna, gdyby można było wykazać, że jest to pewnego rodzaju M-estymator. Niestety nie mogłem wymyślić żadnego, który mógłby dać takie same wyniki.

Oto moja próba. Obliczam logarytmy częstotliwości i próbuję dopasować je do logarytmów teoretycznych prawdopodobieństw podanych przez ten wzór . Ostateczny wynik wydaje się rozsądny. Oto mój kod w 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

Najlepsze dopasowanie kwadratowe to .s=1,47

Maksymalne prawdopodobieństwo w R można wykonać za pomocą mlefunkcji (z stats4pakietu), która pomaga obliczać standardowe błędy (jeśli podano poprawną ujemną funkcję maksymalnego prawdopodobieństwa):

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 

Oto wykres dopasowania do skali log-log (ponownie jak sugerował @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)

Czerwona linia jest sumą dopasowania kwadratów, zielona linia jest dopasowaniem o najwyższym prawdopodobieństwie.

Wykres dziennika pasowań

mpiktas
źródło
1
Jest też pakiet R zipfR cran.r-project.org/web/packages/zipfR/index.html Mimo to nie próbowałem.
onestop
@onstop, dzięki za link. Byłoby miło, gdyby ktoś odpowiedział na to pytanie za pomocą tego pakietu. Moje rozwiązanie zdecydowanie nie ma głębi, choć daje jakąś odpowiedź.
mpiktas
(+1) Jesteś naprawdę imponujący. Tyle dobrych wyników w tak wielu różnych dziedzinach statystyki!
chl
@chl, dzięki! Z pewnością czuję, że nie jestem jedyny z takimi cechami na tej stronie;)
mpiktas
25

Istnieje kilka kwestii przed nami w każdym problemu estymacji:

  1. Oszacuj parametr.

  2. Oceń jakość tego oszacowania.

  3. Przeglądaj dane.

  4. Oceń dopasowanie.

Dla tych, którzy używają metod statystycznych do zrozumienia i komunikacji, pierwsza nigdy nie powinna zostać wykonana bez innych.

ja=1,2),,nja-sss>0

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

ja1n

log(Par(ja))=log(ja-sH.s(n))=-slog(ja)-log(H.s(n)).

faja,ja=1,2),,n

Par(fa1,fa2),,fan)=Par(1)fa1Par(2))fa2)Par(n)fan.

Zatem prawdopodobieństwo dziennika dla danych wynosi

Λ(s)=-sja=1nfajalog(ja)-(ja=1nfaja)log(H.s(n)).

s

s^=1.45041Λ(s^)=-94406,7s^ls=1,463946Λ(s^ls)=-94049,5

s[1,43922,1,46162]

Biorąc pod uwagę naturę prawa Zipf, właściwym sposobem na wykres to dopasowanie jest na wykresie dziennika , gdzie dopasowanie będzie liniowe (z definicji):

wprowadź opis zdjęcia tutaj

Aby ocenić dobroć dopasowania i zbadać dane, spójrz na pozostałości (dane / dopasowanie, osie log-log):

wprowadź opis zdjęcia tutaj

χ2)=656,466


Ponieważ reszty wydają się losowe, w niektórych aplikacjach możemy być zadowoleni z przyjęcia prawa Zipfa (i naszego oszacowania parametru) jako akceptowalnego, choć przybliżonego opisu częstotliwości . Ta analiza pokazuje jednak, że błędem byłoby przypuszczać, że oszacowanie ma jakąkolwiek wartość wyjaśniającą lub predykcyjną dla badanego zestawu danych.

Whuber
źródło
1
@ Whuber, mogę pokornie zasugerować trochę ostrożności w związku z powyższym sformułowaniem. Prawo Zipfa jest zwykle podawane jako wynik częstotliwości względnej. Nie jest to (zwykle uważany) rozkład, z którego pobierana jest próbka. Struktura iid prawdopodobnie nie jest najlepszym pomysłem na te dane. Być może opublikuję więcej na ten temat później.
kardynał
3
@cardinal Czekam na to, co masz do powiedzenia. Jeśli nie masz czasu na dokładną odpowiedź, nawet szkic, który Twoim zdaniem może być „najlepszym pomysłem na te dane”, byłby mile widziany. Mogę zgadnąć, dokąd zmierzasz: dane zostały uszeregowane w rankingu, proces, który tworzy zależności i powinien wymagać ode mnie obrony prawdopodobieństwa uzyskanego bez rozpoznania potencjalnych skutków rankingu. Byłoby miło zobaczyć procedurę szacunkową z uzasadnieniem sygnalizatora. Mam jednak nadzieję, że moją analizę można uratować dzięki rozmiarowi zbioru danych.
whuber
1
@ kardynał, nie rób nam Fermata :) Jeśli masz inny wgląd niż inni, odpowiadaj na nie osobno, nawet jeśli nie będzie to prawidłowa odpowiedź. Na przykład w matematyce takie sytuacje występują dość często.
mpiktas
1
@cardinal łatwo. Na przykład, zbierasz częstotliwości, identyfikujesz i klasyfikujesz dziesięć najwyższych. Podejrzewasz hipotezę prawa Zipfa. Zbierasz nowy zestaw częstotliwości i zgłaszasz je zgodnie z poprzednim rankingiem. To jest sytuacja, do której moja analiza jest idealnie dopasowana, zależnie od tego, czy nowe szeregi zgadzają się ze starymi.
whuber
1
@ whuber, dzięki za cierpliwość. Teraz mam pełną jasność co do twojego rozumowania. Zgodnie z modelem próbkowania, który już w pełni opracowałeś, zgadzam się z twoją analizą. Być może twoje ostatnie stwierdzenie jest wciąż nieco śliskie. Jeśli sortowanie nie wywołuje silnej zależności, wówczas twoja metoda byłaby konserwatywna. Gdyby indukowana zależność była umiarkowanie silna, mogłaby stać się antykonserwatywna. Dziękuję za cierpliwość w obliczu mojej pedanterii.
kardynał
2

s

Jeden z probabilistycznych języków programowania, takich jak PyMC3, sprawia, że ​​oszacowanie jest stosunkowo proste. Inne języki to Stan, który ma wspaniałe funkcje i wspierającą społeczność.

Oto moja implementacja w Pythonie modelu dopasowanego do danych PO (także w 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()

ss

wprowadź opis zdjęcia tutaj

Aby zapewnić podstawową diagnostykę próbkowania, możemy zobaczyć, że próbkowanie było „dobrze mieszane”, ponieważ nie widzimy żadnej struktury w śladzie:

wprowadź opis zdjęcia tutaj

Do uruchomienia kodu potrzebny jest Python z zainstalowanymi pakietami Theano i PyMC3.

Dzięki @ w-huber za wspaniałą odpowiedź i komentarze!

Vladislavs Dovgalecs
źródło
1

Oto moja próba dopasowania danych, oceny i eksploracji wyników za pomocą 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

wprowadź opis zdjęcia tutaj

    Chi-squared test for given probabilities

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

W naszym przypadku zerowymi hipotezami Chi kwadrat jest to, że dane są dystrybuowane zgodnie z prawem zipf, a zatem większe wartości p potwierdzają twierdzenie, że dane są dystrybuowane zgodnie z nim. Zauważ, że nawet bardzo duże wartości p nie są dowodem, a jedynie wskaźnikiem.

Guy s
źródło
0

x=1wx=1^

sUW.S.mi^=H.10-1(1wx=1^)

wx=1^=0,4695599775

sUW.S.mi^=1.4

Ponownie UWSE zapewnia tylko spójne oszacowanie - brak przedziałów ufności, i możemy zobaczyć kompromis w dokładności. Powyższe rozwiązanie mpiktas jest również aplikacją UWSE - chociaż programowanie jest wymagane. Pełne wyjaśnienie estymatora znajduje się na stronie : https://paradsp.wordpress.com/ - na samym dole.

CYP450
źródło
Jak UWSE odnosi się do prawa Zipf?
Michael R. Chernick
UWSE (Unique Weight Space Estimation) wykorzystuje fakt, że najwyższe prawdopodobieństwo / częstotliwość jest unikalne dla różnych wartości parametru s, dla danego N, aby znaleźć s. W odniesieniu do prawa Zipfa mówi nam to, że biorąc pod uwagę liczbę pozycji do uszeregowania, N i najwyższą częstotliwość, istnieje tylko jeden sposób przypisania częstotliwości do pozostałych elementów (2, ..., N), dzięki czemu możemy powiedz „n-ty element jest 1 / n ^ s razy większy niż najczęstszy element, dla niektórych s”. Innymi słowy, biorąc pod uwagę te informacje, istnieje tylko jeden sposób, aby prawo Zipfa utrzymało - oczywiście przy założeniu, że prawo Zipfa faktycznie obowiązuje.
CYP450
0

Moje rozwiązanie stara się być komplementarne do odpowiedzi udzielonych przez mpiktas i whuber wykonujących implementację w Pythonie. Nasze częstotliwości i zakresy x to:

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])

Ponieważ nasza funkcja nie jest zdefiniowana we wszystkich zakresach, musimy sprawdzać, czy normalizujemy się za każdym razem, gdy ją obliczamy. W przypadku dyskretnym prostym przybliżeniem jest podzielenie przez sumę wszystkich y (x). W ten sposób możemy porównać różne parametry.

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)

wprowadź opis zdjęcia tutaj

Wynik daje nam nachylenie 1,450408, jak w poprzednich odpowiedziach.

Ivangtorre
źródło