¿Cómo calcular el coeficiente de la ley de Zipf a partir de un conjunto de frecuencias máximas?

Tengo varias frecuencias de consulta y necesito estimar el coeficiente de la ley de Zipf. Estas son las frecuencias más altas:

26486
12053
5052
3033
2536
2391
1444
1220
1152
1039

Actualización He actualizado el código con el estimador de máxima probabilidad según la sugerencia de @whuber. Reducir al mínimo la suma de los cuadrados de las diferencias entre las probabilidades teóricas del registro y las frecuencias del registro, aunque da una respuesta, sería un procedimiento estadístico si se pudiera demostrar que es algún tipo de estimador M. Lamentablemente, no pude pensar en ninguno que pudiera dar los mismos resultados.

Aquí está mi intento. Calculo los logaritmos de las frecuencias y trato de ajustarlos a los logaritmos de las probabilidades teóricas dadas por esta fórmula . El resultado final parece razonable. Aquí está mi código 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

El mejor ajuste cuadrático entonces es .$s=1.47$

La probabilidad máxima en R se puede realizar con la mlefunción (del stats4paquete), que calcula útilmente los errores estándar (si se proporciona la función de probabilidad máxima negativa correcta):

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 

Aquí está el gráfico del ajuste en la escala log-log (nuevamente como sugirió @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)

La línea roja es la suma de los cuadrados ajustados, la línea verde es el ajuste de máxima probabilidad.

Hay varios problemas ante nosotros en cualquier problema de estimación:

1. Estime el parámetro.

2. Evaluar la calidad de esa estimación.

3. Explore los datos.

4. Evaluar el ajuste.

Para aquellos que usarían métodos estadísticos para la comprensión y la comunicación, el primero nunca debe hacerse sin los demás.

$i=1,2,\ldots,n$$i^{-s}$$s$$s\gt 0$

$i$$1$$n$

$f_i, i=1,2,\ldots, n$

Por lo tanto, la probabilidad de registro para los datos es

$s$ .

$\hat{s} = 1.45041$$\Lambda(\hat{s}) = -94046.7$$\hat{s}_{ls} = 1.463946$$\Lambda(\hat{s}_{ls}) = -94049.5$

ML también estimará $s$$[1.43922, 1.46162]$ (si hice los cálculos correctamente :-).

Dada la naturaleza de la ley de Zipf, la forma correcta de graficar este ajuste es en un gráfico log-log , donde el ajuste será lineal (por definición):

Para evaluar la bondad del ajuste y explorar los datos, observe los residuos (datos / ajuste, ejes log-log nuevamente):

Esto no es demasiado bueno: aunque no hay una correlación serial evidente o heteroscedasticidad en los residuos, generalmente son alrededor del 10% (lejos de 1.0). Con frecuencias de miles, no esperaríamos desviaciones en más de un pequeño porcentaje. los$\chi^2 = 656.476$

Debido a que los residuos parecen aleatorios, en algunas aplicaciones podríamos contentarnos con aceptar la Ley de Zipf (y nuestra estimación del parámetro) como una descripción aceptable aunque aproximada de las frecuencias . Sin embargo, este análisis muestra que sería un error suponer que esta estimación tiene algún valor explicativo o predictivo para el conjunto de datos examinado aquí.

Las estimaciones de máxima verosimilitud son solo estimaciones puntuales del parámetro $s$ . Se necesita un esfuerzo adicional para encontrar también el intervalo de confianza de la estimación. El problema es que este intervalo no es probabilístico. No se puede decir "el valor del parámetro s = ... tiene una probabilidad del 95% en el rango [...]".

Uno de los lenguajes de programación probabilísticos como PyMC3 hace que esta estimación sea relativamente sencilla. Otros idiomas incluyen Stan, que tiene excelentes características y una comunidad de apoyo.

Aquí está mi implementación de Python del modelo ajustado en los datos de los OP (también en 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()

$s$$s$

Para proporcionar algunos diagnósticos básicos de muestreo, podemos ver que el muestreo estaba "mezclando bien" ya que no vemos ninguna estructura en la traza:

Para ejecutar el código, se necesita Python con los paquetes Theano y PyMC3 instalados.

¡Gracias a @ w-huber por su excelente respuesta y comentarios!

Aquí está mi intento de ajustar los datos, evaluar y explorar los resultados usando 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

En nuestro caso, las hipótesis nulas de Chi cuadrado son que los datos se distribuyen de acuerdo con la ley de zipf, por lo tanto, los valores p más grandes respaldan la afirmación de que los datos se distribuyen de acuerdo con ellos. Tenga en cuenta que incluso los valores p muy grandes no son una prueba, solo un indicador.

$x = 1$$\hat{w_{x=1}}$

$\hat{w_{x=1}} = 0.4695599775$

Nuevamente, el UWSE solo proporciona una estimación consistente, sin intervalos de confianza, y podemos ver cierta compensación en la precisión. La solución anterior de mpiktas también es una aplicación de UWSE, aunque se requiere programación. Para obtener una explicación completa del estimador, consulte: https://paradsp.wordpress.com/ - hasta el final.

Mi solución intenta ser complementaria a las respuestas proporcionadas por mpiktas y whuber haciendo una implementación en Python. Nuestras frecuencias y rangos x son:

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 nuestra función no está definida en todo el rango, debemos verificar que estamos normalizando cada vez que la calculamos. En el caso discreto, una aproximación simple es dividir por la suma de todos y (x). De esta 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)

El resultado nos da una pendiente de 1.450408 como en las respuestas anteriores.