Como calcular o coeficiente da lei de Zipf a partir de um conjunto de frequências mais altas?

25

Tenho várias frequências de consulta e preciso estimar o coeficiente da lei de Zipf. Estas são as principais frequências:

26486
12053
5052
3033
2536
2391
1444
1220
1152
1039
Diegolo
fonte
de acordo com a página da wikipedia, a lei de Zipf tem dois parâmetros. Número dos elementos e o expoente. O que é no seu caso, 10? E as frequências podem ser calculadas dividindo os valores fornecidos pela soma de todos os valores fornecidos? s NNsN
Mvctas
deixe dez e as frequências podem ser calculadas dividindo os valores fornecidos pela soma de todos os valores fornecidos. como posso estimar?
Diegolo

Respostas:

22

Atualização Atualizei o código com o estimador de máxima verossimilhança, conforme a sugestão @whuber. Minimizar a soma dos quadrados das diferenças entre as probabilidades teóricas do log e as frequências do log, embora dê uma resposta, seria um procedimento estatístico se fosse possível demonstrar que é algum tipo de estimador-M. Infelizmente, não consegui pensar em nenhum que pudesse dar os mesmos resultados.

Aqui está a minha tentativa. Calculo os logaritmos das frequências e tento ajustá-los aos logaritmos das probabilidades teóricas fornecidas por esta fórmula . O resultado final parece razoável. Aqui está o meu código em 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

O melhor ajuste quadrático então é .s=1.47

A probabilidade máxima em R pode ser executada com a mlefunção (do stats4pacote), que calcula os erros padrão de maneira útil (se a função de probabilidade máxima negativa correta for fornecida):

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 

Aqui está o gráfico do ajuste na escala log-log (novamente como @whuber sugeriu):

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)

A linha vermelha é a soma dos quadrados, a linha verde é a probabilidade máxima.

Gráfico de log-log de ajustes

mpiktas
fonte
11
Há também um pacote R zipfR cran.r-project.org/web/packages/zipfR/index.html Eu ainda não tentei.
onestop
@ onestop, obrigado pelo link. Seria bom se alguém respondesse a essa pergunta usando este pacote. Minha solução definitivamente não tem profundidade, embora dê algum tipo de resposta.
Mpiktas
(+1) Você é realmente impressionante. Tantas boas contribuições em tantos campos estatísticos diferentes!
chl
@chl, obrigado! Eu certamente sinto que eu não sou o único com tais características deste site embora;)
mpiktas
25

Existem vários problemas diante de nós em qualquer problema de estimativa:

  1. Estime o parâmetro.

  2. Avalie a qualidade dessa estimativa.

  3. Explore os dados.

  4. Avalie o ajuste.

Para aqueles que usariam métodos estatísticos para compreensão e comunicação, o primeiro nunca deveria ser feito sem os outros.

Para a estimativa , é conveniente usar a máxima verossimilhança (ML). As frequências são tão grandes que podemos esperar que as conhecidas propriedades assintóticas se mantenham. ML usa a distribuição de probabilidade assumida dos dados. A lei de Zipf supõe que as probabilidades de sejam proporcionais a para algumas potências constantes (geralmente ). Como essas probabilidades devem somar à unidade, a constante de proporcionalidade é o recíproco da somai=1,2,,nisss>0

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

Consequentemente, o logaritmo da probabilidade de qualquer resultado entre e éi1n

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

Para dados independentes resumidos por suas frequências , a probabilidade é o produto das probabilidades individuais,fi,i=1,2,,n

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

Assim, a probabilidade do log para os dados é

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

Considerando os dados como fixos e expressando-os explicitamente como uma função de , torna-se o log Probabilidade .s

Minimização numérica do log Probabilidade com os dados fornecidos na pergunta produz e . Isso é significativamente melhor (mas apenas pouco) do que a solução dos mínimos quadrados (com base nas frequências de log) de com . (A otimização pode ser feita com uma pequena alteração no código R claro e elegante fornecido por mpiktas.)s^=1.45041Λ(s^)=94046.7s^ls=1.463946Λ(s^ls)=94049.5

O ML também estimará os limites de confiança para da maneira usual. A aproximação do qui-quadrado fornece (se eu fiz os cálculos corretamente :-).[ 1.43922 , 1,46162 ]s[1.43922,1.46162]

Dada a natureza da lei de Zipf, a maneira correta de representar graficamente esse ajuste é em um gráfico de log-log , em que o ajuste será linear (por definição):

insira a descrição da imagem aqui

Para avaliar a qualidade do ajuste e explorar os dados, observe os resíduos (dados / ajuste, eixos de log-log novamente):

insira a descrição da imagem aqui

Isso não é muito grande: embora não exista correlação serial evidente ou heterocedasticidade nos resíduos, eles geralmente estão em torno de 10% (longe de 1,0). Com frequências na casa dos milhares, não esperaríamos desvios em mais do que alguns por cento. A qualidade do ajuste é prontamente testada com o qui quadrado . Obtemos com = 9 graus de liberdade; isso é uma evidência altamente significativa de desvios da lei de Zipf .χ2=656.476


Como os resíduos parecem aleatórios, em algumas aplicações, podemos aceitar aceitar a Lei de Zipf (e nossa estimativa do parâmetro) como uma descrição aceitável, embora grosseira, das frequências . Essa análise mostra, no entanto, que seria um erro supor que essa estimativa tenha algum valor explicativo ou preditivo para o conjunto de dados examinado aqui.

whuber
fonte
11
@ Whuber, posso humildemente sugerir um pouco de cautela com a formulação dada acima. A lei de Zipf é geralmente declarada como um resultado de frequência relativa. Não é (normalmente considerada) a distribuição a partir da qual uma amostra iid é retirada. Uma estrutura iid provavelmente não é a melhor ideia para esses dados. Talvez eu poste mais sobre isso mais tarde.
cardinal
3
@ cardinal Estou ansioso pelo que você tem a dizer. Se você não tiver tempo para uma resposta completa, mesmo um esboço do que você acha que pode ser a "melhor idéia para esses dados" será bem-vindo. Posso adivinhar aonde você está indo com isso: os dados foram classificados, um processo que cria dependências e deve exigir que eu defenda uma probabilidade derivada sem reconhecer os efeitos potenciais do ranking. Seria bom ver um procedimento de estimativa com justificativa mais sólida. Estou esperançoso, no entanto, que minha análise possa ser resgatada pelo tamanho do conjunto de dados.
whuber
11
@ cardinal, não faça Fermat conosco :) Se você tiver uma visão diferente da de outros respondentes, sinta-se à vontade para expressá-la na resposta separada, mesmo que ela não constitua uma resposta válida por si só. Em math.SE, por exemplo, tais situações surgem com bastante frequência.
Mpiktas
11
@ cardinal Facilmente. Por exemplo, você coleta frequências e identifica e classifica as dez mais altas. Você supõe a Lei de Zipf. Você coleta um novo conjunto de frequências e as informa de acordo com a classificação anterior . Essa é a situação da segunda, à qual minha análise é perfeitamente adequada, dependendo das novas fileiras concordando com as antigas.
whuber
11
@ whuber, obrigado por sua paciência. Agora estou totalmente claro sobre sua linha de raciocínio. Sob o modelo de amostragem que você desenvolveu agora, concordo com sua análise. Talvez sua última declaração ainda seja um pouco escorregadia. Se a classificação não induzir forte dependência, seu método seria conservador. Se a dependência induzida for moderadamente forte, ela poderá se tornar anticonservadora. Agradecemos sua paciência diante do meu pedantismo.
cardeal
2

As estimativas de máxima verossimilhança são apenas estimativas pontuais dos parâmetros . É necessário um esforço extra para encontrar também o intervalo de confiança da estimativa. O problema é que esse intervalo não é probabilístico. Não se pode dizer "o valor do parâmetro s = ... está com probabilidade de 95% na faixa [...]".s

Uma das linguagens de programação probabilística, como PyMC3, torna essa estimativa relativamente direta. Outros idiomas incluem Stan, que possui ótimos recursos e comunidade de suporte.

Aqui está minha implementação em Python do modelo ajustado nos dados dos OPs (também no 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()

Aqui estão as estimativas dos parâmetros na forma de distribuição. Observe como compacta é a estimativa! Com probabilidade de 95%, o valor verdadeiro do parâmetro está na faixa [1.439,1.461]; a média é de cerca de 1,45, o que está muito próximo das estimativas do MLE.sss

insira a descrição da imagem aqui

Para fornecer alguns diagnósticos básicos de amostragem, podemos ver que a amostragem estava "se misturando bem", pois não vemos nenhuma estrutura no rastreamento:

insira a descrição da imagem aqui

Para executar o código, é necessário o Python com os pacotes Theano e PyMC3 instalados.

Obrigado a @ w-huber por sua ótima resposta e comentários!

Vladislavs Dovgalecs
fonte
1

Aqui está minha tentativa de ajustar os dados, avaliar e explorar os resultados usando o 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

insira a descrição da imagem aqui

    Chi-squared test for given probabilities

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

No nosso caso, as hipóteses nulas de Chi square são que os dados são distribuídos de acordo com a lei do zipf, portanto, valores p maiores sustentam a alegação de que os dados são distribuídos de acordo com ele. Observe que mesmo valores p muito grandes não são uma prova, apenas um indicador.

Rapazes
fonte
0

Apenas por diversão, essa é outra instância em que o UWSE pode fornecer uma solução de formulário fechado usando apenas a maior frequência - embora a um custo de precisão. A probabilidade em é única entre os valores dos parâmetros. Se denota a frequência relativa correspondente, então,x=1wx=1^

sUWSE^=H101(1wx=1^)

Nesse caso, como , obtemos:wx=1^=0.4695599775

sUWSE^=1.4

Novamente, o UWSE fornece apenas uma estimativa consistente - sem intervalos de confiança, e podemos ver algumas compensações na precisão. A solução do mpiktas acima também é uma aplicação do UWSE - embora seja necessária programação. Para uma explicação completa do estimador, consulte: https://paradsp.wordpress.com/ - todo o caminho na parte inferior.

CYP450
fonte
Como o UWSE se relaciona com a lei de Zipf?
Michael R. Chernick
O UWSE (Estimativa de espaço de peso exclusivo) usa o fato de que a probabilidade / frequência máxima é única entre diferentes valores do parâmetro s, para um dado N, para encontrar s. Com relação à lei de Zipf, isso nos diz que, dado um número de itens para classificar, N e uma frequência máxima, há apenas uma maneira de atribuir frequências aos itens restantes (2, ..., N), para que possamos diga "o enésimo item é 1 / n ^ s vezes maior que o item mais frequente, para alguns s". Em outras palavras, dada essa informação, existe apenas uma maneira de a lei de Zipf se manter - é claro, supondo que a lei de Zipf se mantenha.
precisa
0

Minha solução tenta complementar as respostas fornecidas por mpiktas e whuber fazendo uma implementação em Python. Nossas frequências e faixas x são:

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

Como nossa função não está definida em todos os intervalos, precisamos verificar se estamos normalizando cada vez que a computamos. No caso discreto, uma aproximação simples é dividir pela soma de todos os y (x). Desta forma, podemos comparar diferentes parâmetros.

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)

insira a descrição da imagem aqui

O resultado nos dá uma inclinação de 1.450408, como nas respostas anteriores.

ivangtorre
fonte