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
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
Réponses:
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.
Le meilleur ajustement quadratique est alors .s = 1,47
La probabilité maximale dans R peut être effectuée avec la
mle
fonction (dustats4
package), qui calcule utilement les erreurs standard (si la fonction de probabilité maximale négative correcte est fournie):Voici le graphique de l'ajustement dans l'échelle log-log (encore une fois comme @whuber l'a suggéré):
La ligne rouge correspond à la somme des carrés, la ligne verte correspond à l'ajustement de probabilité maximale.
la source
Il y a plusieurs problèmes devant nous dans tout problème d'estimation:
Estimez le paramètre.
Évaluez la qualité de cette estimation.
Explorez les données.
É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.i = 1 , 2 , … , n je- s s s > 0
Ainsi, la probabilité logarithmique des données est
É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):
Pour évaluer la qualité de l'ajustement et explorer les données, regardez les résidus (données / ajustement, axes log-log à nouveau):
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.
la source
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 ):
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:
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!
la source
Voici ma tentative d'adapter les données, d'évaluer et d'explorer les résultats à l'aide de VGAM:
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.
la source
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 = 1 est unique parmi les valeurs des paramètres. Siwx = 1^ désigne alors la fréquence relative correspondante,
Dans ce cas, puisquewx = 1^= 0,4695599775 , on a:
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.
la source
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:
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.
Le résultat nous donne une pente de 1,450408 comme dans les réponses précédentes.
la source