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

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

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.

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 soma$i=1,2,\ldots,n$$i^{-s}$$s$$s\gt 0$

Consequentemente, o logaritmo da probabilidade de qualquer resultado entre e é$i$$1$$n$

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

Assim, a probabilidade do log para os dados é

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.)$\hat{s} = 1.45041$$\Lambda(\hat{s}) = -94046.7$$\hat{s}_{ls} = 1.463946$$\Lambda(\hat{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):

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

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 .$\chi^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.

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.$s$$s$

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:

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!

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

    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.

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 = 1$$\hat{w_{x=1}}$

Nesse caso, como , obtemos:$\hat{w_{x=1}} = 0.4695599775$

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.

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)

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