Comment calculer le coefficient de loi de Zipf à partir d'un ensemble de fréquences maximales?

25

J'ai plusieurs fréquences de requête et j'ai besoin d'estimer le coefficient de la loi de Zipf. Ce sont les fréquences les plus élevées:

26486
12053
5052
3033
2536
2391
1444
1220
1152
1039
Diegolo
la source
selon la page wikipedia, la loi de Zipf a deux paramètres. Nombre des éléments N et s l'exposant. Quel est N dans votre cas, 10? Et les fréquences peuvent être calculées en divisant vos valeurs fournies par la somme de toutes les valeurs fournies?
mpiktas
laissez-le dix, et les fréquences peuvent être calculées en divisant vos valeurs fournies par la somme de toutes les valeurs fournies .. comment puis-je estimer?
Diegolo

Réponses:

22

Mettre à jour J'ai mis à jour le code avec l'estimateur du maximum de vraisemblance selon la suggestion de @whuber. Minimiser la somme des carrés des différences entre les probabilités théoriques logarithmiques et les fréquences logarithmiques donne une réponse serait une procédure statistique s'il pouvait être démontré qu'il s'agit d'une sorte d'estimateur M. Malheureusement, je n'ai pu penser à aucun qui pourrait donner les mêmes résultats.

Voici ma tentative. Je calcule les logarithmes des fréquences et essaie de les adapter aux logarithmes des probabilités théoriques données par cette formule . Le résultat final semble raisonnable. Voici mon code en 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

Le meilleur ajustement quadratique est alors .s=1,47

La probabilité maximale dans R peut être effectuée avec la mlefonction (du stats4package), qui calcule utilement les erreurs standard (si la fonction de probabilité maximale négative correcte est fournie):

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 

Voici le graphique de l'ajustement dans l'échelle log-log (encore une fois comme @whuber l'a suggéré):

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)

La ligne rouge correspond à la somme des carrés, la ligne verte correspond à l'ajustement de probabilité maximale.

Log-log graphique d'ajustements

mpiktas
la source
1
Il y a aussi un package R zipfR cran.r-project.org/web/packages/zipfR/index.html je ne l'ai pas essayé cependant.
onestop
@onestop, merci pour le lien. Ce serait bien si quelqu'un répondait à cette question en utilisant ce package. Ma solution manque définitivement de profondeur, bien qu'elle donne une sorte de réponse.
mpiktas
(+1) Vous êtes vraiment impressionnant. Tant de bonnes contributions dans tant de domaines statistiques différents!
chl
@chl, merci! Je pense certainement que je ne suis pas le seul à avoir de telles caractéristiques sur ce site;)
mpiktas
25

Il y a plusieurs problèmes devant nous dans tout problème d'estimation:

  1. Estimez le paramètre.

  2. Évaluez la qualité de cette estimation.

  3. Explorez les données.

  4. Évaluez l'ajustement.

Pour ceux qui utiliseraient des méthodes statistiques pour comprendre et communiquer, la première ne devrait jamais se faire sans les autres.

Pour l' estimation, il est pratique d'utiliser la vraisemblance maximale (ML). Les fréquences sont si grandes que nous pouvons nous attendre à ce que les propriétés asymptotiques bien connues se maintiennent. ML utilise la distribution de probabilité supposée des données.je=1,2,,nje-sss>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.

Ainsi, la probabilité logarithmique des données est

Λ(s)=-sje=1nFjebûche(je)-(je=1nFje)bûche(Hs(n)).

s

s^=1.45041Λ(s^)=-94046.7s^ls=1,463946Λ(s^ls)=-94049.5

s[1.43922,1,46162]

Étant donné la nature de la loi de Zipf, la bonne façon de représenter graphiquement cet ajustement est sur un tracé log-log , où l'ajustement sera linéaire (par définition):

entrez la description de l'image ici

Pour évaluer la qualité de l'ajustement et explorer les données, regardez les résidus (données / ajustement, axes log-log à nouveau):

entrez la description de l'image ici

χ2=656.476


Parce que les résidus semblent aléatoires, dans certaines applications, nous pourrions nous contenter d'accepter la loi de Zipf (et notre estimation du paramètre) comme une description acceptable mais grossière des fréquences . Cette analyse montre cependant que ce serait une erreur de supposer que cette estimation a une valeur explicative ou prédictive pour l'ensemble de données examiné ici.

whuber
la source
1
@whuber, je pourrais humblement suggérer un peu de prudence avec la formulation donnée ci-dessus. La loi de Zipf est généralement indiquée comme un résultat de fréquence relative. Ce n'est pas (normalement considéré) la distribution à partir de laquelle un échantillon iid est tiré. Un cadre iid n'est probablement pas la meilleure idée pour ces données. Peut-être que je posterai plus à ce sujet plus tard.
Cardinal
3
@cardinal J'attends avec impatience ce que vous avez à dire. Si vous n'avez pas le temps pour une réponse approfondie, même une esquisse de ce que vous pensez être la "meilleure idée pour ces données" serait la bienvenue. Je peux deviner où vous allez avec ceci: les données ont été classées, un processus qui crée des dépendances et devrait m'obliger à défendre une probabilité dérivée sans reconnaître les effets potentiels du classement. Ce serait bien de voir une procédure d'estimation avec une justification plus solide. J'espère cependant que mon analyse pourra être sauvée par la taille même de l'ensemble de données.
whuber
1
@cardinal, ne faites pas Fermat sur nous :) Si vous avez un aperçu différent de celui des autres répondants, n'hésitez pas à l'exprimer dans la réponse séparée, même si cela ne constituera pas une réponse valide en soi. En math.SE par exemple, de telles situations se produisent assez fréquemment.
mpiktas
1
@cardinal facilement. Par exemple, vous collectez des fréquences et identifiez et classez les dix plus hautes. Vous faites l'hypothèse de la loi de Zipf. Vous collectez un nouvel ensemble de fréquences et les rapportez selon le classement précédent . C'est la situation iid, à laquelle mon analyse est parfaitement adaptée, sous réserve que les nouveaux rangs soient d'accord avec les anciens.
whuber
1
@whuber, merci pour votre patience. Maintenant, je comprends parfaitement votre raisonnement. Dans le modèle d'échantillonnage que vous avez maintenant entièrement étoffé, je suis d'accord avec votre analyse. Votre toute dernière déclaration est peut-être encore un peu glissante. Si le tri n'induit pas une forte dépendance, votre méthode serait conservatrice. Si la dépendance induite était modérément forte, elle pourrait devenir anticonservatrice. Merci pour votre patience face à mon pédantisme.
Cardinal
2

s

L'un des langages de programmation probabilistes comme PyMC3 rend cette estimation relativement simple. D'autres langues incluent Stan qui a de grandes fonctionnalités et une communauté de soutien.

Voici mon implémentation Python du modèle monté sur les données OPs (également sur 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()

ssous forme de distribution. Remarquez à quel point l'estimation est compacte! Avec une probabilité de 95%, la vraie valeur du paramètresest dans la plage [1.439,1.461]; la moyenne est d'environ 1,45, ce qui est très proche des estimations MLE.

entrez la description de l'image ici

Pour fournir des diagnostics d'échantillonnage de base, nous pouvons voir que l'échantillonnage "se mélangeait bien" car nous ne voyons aucune structure dans la trace:

entrez la description de l'image ici

Pour exécuter le code, il faut Python avec les packages Theano et PyMC3 installés.

Merci à @ w-huber pour sa grande réponse et ses commentaires!

Vladislavs Dovgalecs
la source
1

Voici ma tentative d'adapter les données, d'évaluer et d'explorer les résultats à l'aide de 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

entrez la description de l'image ici

    Chi-squared test for given probabilities

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

Dans notre cas, l'hypothèse nulle du chi carré est que les données sont distribuées selon la loi de zipf, donc des valeurs de p plus grandes soutiennent l'affirmation selon laquelle les données sont distribuées selon elle. Notez que même de très grandes valeurs de p ne sont pas une preuve, juste un indicateur.

Les gars
la source
0

Juste pour le plaisir, c'est un autre exemple où l'UWSE peut fournir une solution sous forme fermée utilisant uniquement la fréquence la plus élevée - mais à un coût de précision. La probabilité surX=1est unique parmi les valeurs des paramètres. SiwX=1^ désigne alors la fréquence relative correspondante,

sUWSE^=Hdix-1(1wX=1^)

Dans ce cas, puisque wX=1^=0,4695599775, on a:

sUWSE^=1.4

Encore une fois, l'UWSE ne fournit qu'une estimation cohérente - aucun intervalle de confiance, et nous pouvons voir un certain compromis dans la précision. La solution de mpiktas ci-dessus est également une application de l'UWSE - bien que la programmation soit requise. Pour une explication complète de l'estimateur, voir: https://paradsp.wordpress.com/ - tout en bas.

CYP450
la source
Quel est le lien entre l'UWSE et la loi de Zipf?
Michael R. Chernick
L'UWSE (Unique Weight Space Estimation) utilise le fait que la probabilité / fréquence la plus élevée est unique à travers différentes valeurs du paramètre s, pour un N donné, pour trouver s. En ce qui concerne la loi de Zipf, cela nous dit que, étant donné un certain nombre d'éléments à classer, N et une fréquence maximale, il n'y a qu'une seule façon d'attribuer des fréquences aux éléments restants (2, ..., N) de sorte que nous puissions dire "le nième élément est 1 / n ^ fois plus grand que l'élément le plus fréquent, pour certains s". En d'autres termes, compte tenu de ces informations, il n'y a qu'une seule façon pour la loi de Zipf de tenir - bien sûr, en supposant que la loi de Zipf est effectivement valable.
CYP450
0

Ma solution essaie d'être complémentaire aux réponses fournies par mpiktas et whuber faisant une implémentation en Python. Nos fréquences et gammes x sont:

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

Comme notre fonction n'est pas définie dans toutes les plages, nous devons vérifier que nous normalisons chaque fois que nous la calculons. Dans le cas discret, une approximation simple consiste à diviser par la somme de tous les y (x). De cette façon, nous pouvons comparer différents paramètres.

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)

entrez la description de l'image ici

Le résultat nous donne une pente de 1,450408 comme dans les réponses précédentes.

ivangtorre
la source