Как рассчитать коэффициент закона Ципфа из набора верхних частот?

25

У меня есть несколько частот запросов, и мне нужно оценить коэффициент закона Ципфа. Это верхние частоты:

26486
12053
5052
3033
2536
2391
1444
1220
1152
1039
Diegolo
источник
согласно странице википедии закон Ципфа имеет два параметра. Количество элементов и степени. Что такое в вашем случае, 10? А частоты можно рассчитать путем деления предоставленных вами значений на сумму всех предоставленных значений? сек NNsN
mpiktas
пусть это десять, и частоты могут быть рассчитаны путем деления предоставленных вами значений на сумму всех предоставленных значений. Как я могу оценить?
Diegolo

Ответы:

22

Обновление Я обновил код с оценкой максимального правдоподобия согласно предложению @whuber. Минимизация суммы квадратов различий между логарифмическими теоретическими вероятностями и логарифмическими частотами, хотя дает ответ, была бы статистической процедурой, если бы можно было показать, что это своего рода М-оценка. К сожалению, я не мог придумать ни одного, который мог бы дать такие же результаты.

Вот моя попытка. Я вычисляю логарифмы частот и пытаюсь привести их в соответствие с логарифмами теоретических вероятностей, которые дает эта формула . Конечный результат кажется разумным. Вот мой код в 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

Тогда наилучшее квадратичное соответствие будет .sзнак равно1,47

Максимальное правдоподобие в R можно выполнить с помощью mleфункции (из stats4пакета), которая услужливо рассчитывает стандартные ошибки (если указана правильная функция отрицательного максимального правдоподобия):

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 

Вот график соответствия в масштабе log-log (опять же, как предложил @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)

Красная линия - это сумма подходящих квадратов, зеленая линия - соответствие с максимальным правдоподобием.

Лог-график граф посадки

mpiktas
источник
1
Есть также пакет R zipfR cran.r-project.org/web/packages/zipfR/index.html, хотя я еще не пробовал.
OneStop
@onestop, спасибо за ссылку. Было бы хорошо, если бы кто-то ответил на этот вопрос, используя этот пакет. Моему решению определенно не хватает глубины, хотя оно дает какой-то ответ.
mpiktas
(+1) Ты действительно впечатляет. Так много хороших вкладов в самых разных областях статистики!
хл
@ CHL, спасибо! Я, конечно, чувствую, что я не единственный с такими характеристиками на этом сайте;)
mpiktas
25

Перед любой проблемой оценки стоит несколько вопросов :

  1. Оцените параметр.

  2. Оцените качество этой оценки.

  3. Изучите данные.

  4. Оцените подгонку.

Для тех, кто использует статистические методы для понимания и общения, первое никогда не должно делаться без других.

язнак равно1,2,...,Nя-sss>0

ЧАСs(N)знак равно11s+12s++1Ns,

я1N

журнал(Pr(я))знак равножурнал(я-sЧАСs(N))знак равно-sжурнал(я)-журнал(ЧАСs(N)),

ея,язнак равно1,2,...,N

Pr(е1,е2,...,еN)знак равноPr(1)е1Pr(2)е2Pr(N)еN,

Таким образом, логарифмическая вероятность для данных

Λ(s)знак равно-sΣязнак равно1Nеяжурнал(я)-(Σязнак равно1Nея)журнал(ЧАСs(N)),

s

s^знак равно1,45041Λ(s^)знак равно-94046,7s^Lsзнак равно1.463946Λ(s^Ls)знак равно-94049,5

s[1,43922,1,46162]

Учитывая природу закона Ципфа, правильный способ составить график этого подбора - на графике log-log , где подгонка будет линейной (по определению):

введите описание изображения здесь

Чтобы оценить правильность подбора и изучить данные, посмотрите на остатки (данные / подбор, снова ось log-log):

введите описание изображения здесь

χ2знак равно656,476


Поскольку остатки кажутся случайными, в некоторых приложениях мы могли бы принять закон Ципфа (и нашу оценку параметра) в качестве приемлемого, хотя и грубого описания частот . Этот анализ показывает, однако, что было бы ошибкой предполагать, что эта оценка имеет какое-либо объяснительное или прогнозное значение для набора данных, рассматриваемого здесь.

Whuber
источник
1
@whuber, я мог бы смиренно предложить немного осторожности с формулировкой, приведенной выше. Закон Ципфа обычно указывается как результат с относительной частотой. Это не (обычно считается) распределение, из которого берется образец iid. Фреймворк iid, вероятно, не лучшая идея для этих данных. Возможно, я опубликую больше об этом позже.
кардинал
3
@cardinal Я с нетерпением жду, что вы скажете. Если у вас нет времени для тщательного ответа, даже набросок того, что, по вашему мнению, может быть «лучшей идеей для этих данных», будет очень кстати. Я могу догадаться, куда вы идете с этим: данные были ранжированы, процесс, который создает зависимости и должен требовать от меня защиты вероятности, полученной без признания потенциальных эффектов ранжирования. Было бы неплохо увидеть процедуру оценки с более обоснованным обоснованием. Тем не менее, я надеюсь, что мой анализ может быть спасен огромным размером набора данных.
whuber
1
@cardinal, не делайте Ферма на нас :) Если у вас есть какое-то иное понимание, чем у других отвечающих, не стесняйтесь выразить его в отдельном ответе, даже если он не будет являться действительным ответом как таковым. В math.SE , например , возникают такие ситуации довольно часто.
mpiktas
1
@ Cardinal Легко. Например, вы собираете частоты, идентифицируете и оцениваете десятку самых высоких. Вы выдвигаете гипотезу о законе Ципфа. Вы собираете новый набор частот и сообщаете о них в соответствии с предыдущим рейтингом. Это та ситуация, в которой мой анализ идеально подходит, в зависимости от того, согласны ли новые ряды со старыми.
whuber
1
@ whuber, спасибо за ваше терпение. Теперь я полностью понимаю вашу аргументацию. В соответствии с моделью выборки, которую вы сейчас полностью определили, я согласен с вашим анализом. Возможно, ваше последнее заявление все еще немного скользко. Если сортировка не вызывает сильной зависимости, тогда ваш метод будет консервативным. Если бы индуцированная зависимость была умеренно сильной, она могла бы стать антиконсервативной. Спасибо за терпение перед лицом моей педантичности.
кардинал
2

s

Один из вероятностных языков программирования, такой как PyMC3, делает эту оценку относительно простой. Другие языки включают Stan, у которого есть отличные возможности и поддерживающее сообщество.

Вот моя реализация Python модели, адаптированной к данным OP (также на 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()

ss

введите описание изображения здесь

Чтобы обеспечить некоторую базовую диагностику выборки, мы можем видеть, что выборка «хорошо перемешивалась», так как мы не видим никакой структуры в трассе:

введите описание изображения здесь

Для запуска кода необходим Python с установленными пакетами Theano и PyMC3.

Спасибо @ w-huber за отличный ответ и комментарии!

Владислав Довгальец
источник
1

Вот моя попытка привести данные в соответствие, оценить и изучить результаты с помощью 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

В нашем случае нулевая гипотеза Чи-квадрата состоит в том, что данные распределяются по закону Зипфа, поэтому большие значения р подтверждают утверждение, что данные распределяются в соответствии с ним. Обратите внимание, что даже очень большие значения p не являются доказательством, это просто показатель.

Парень с
источник
0

Иксзнак равно1весИксзнак равно1^

sUWSЕ^знак равноЧАС10-1(1весИксзнак равно1^)

весИксзнак равно1^знак равно0,4695599775

sUWSЕ^знак равно1.4

Опять же, UWSE предоставляет только непротиворечивую оценку - без доверительных интервалов, и мы можем увидеть некоторый компромисс в точности. Приведенное выше решение mpiktas также является приложением UWSE - хотя программирование необходимо. Для полного объяснения оценки см .: https://paradsp.wordpress.com/ - полностью внизу.

CYP450
источник
Как UWSE относится к закону Ципфа?
Майкл Р. Черник
UWSE (Уникальная оценка пространства весов) использует тот факт, что самая высокая вероятность / частота уникальна для разных значений параметра s для данного N, чтобы найти s. Что касается закона Ципфа, это говорит нам о том, что, учитывая количество элементов для ранжирования, N и самую верхнюю частоту, существует только один способ присвоения частот оставшимся элементам (2, ..., N), так что мы можем скажем "n-й элемент в 1 / n ^ s раз больше, чем самый частый элемент, для некоторых s". Другими словами, учитывая эту информацию, закон Ципфа может существовать только одним способом - конечно, предполагая, что закон Ципфа действительно действует.
CYP450
0

Мое решение пытается дополнить ответы, предоставленные mpiktas и whuber, которые осуществляют реализацию на Python. Наши частоты и диапазоны х:

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

Поскольку наша функция не определена во всем диапазоне, нам нужно проверять, что мы нормализуем каждый раз, когда вычисляем ее. В дискретном случае простым приближением является деление на сумму всех y (x). Таким образом, мы можем сравнить различные параметры.

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)

введите описание изображения здесь

Результат дает нам наклон 1.450408, как и в предыдущих ответах.

ivangtorre
источник