Как оценить параметры для усеченного распределения Zipf из выборки данных?

У меня проблема с параметром оценки для Zipf. Моя ситуация следующая:

У меня есть набор образцов (измеренный из эксперимента, который генерирует вызовы, которые должны следовать распределению Zipf). Я должен продемонстрировать, что этот генератор действительно генерирует вызовы с распределением zipf. Я уже читал эти вопросы и ответы. Как рассчитать коэффициент закона Ципфа из набора верхних частот? но я достигаю плохих результатов, потому что я использую усеченный дистрибутив. Например, если я устанавливаю значение «s» на «0,9» для процесса генерации, если я пытаюсь оценить значение «s», как написано в сообщенных вопросах и ответах, я получаю «s», равный 0,2 ca. Я думаю, это связано с тем, что я использую дистрибутив TRUNCATED (мне нужно ограничить zipf точкой усечения, он усекается вправо).

Как я могу оценить параметры с усеченным распределением zipf?

distributions estimation pareto-distribution zipf Maurizio
источник

чтобы быть ясно, что именно вы правы усечения? Распределение значений или сам график Zipf? Вы знаете точку усечения? Является ли усечение артефактом данных или артефактом обработки данных (например, какое-то решение, принятое вами или экспериментатором)? Любые дополнительные детали будут полезны.

кардинал

@cardinal. (часть 1/2) Спасибо кардинал. Я дам более подробную информацию: у меня есть генератор VoIP, который генерирует вызовы после Zipf (и другого распределения) для объема на одного абонента. Я должен убедиться, что этот генератор действительно следует этим распределениям. Для Zipf Distribution мне нужно определить точку усечения (следовательно, она известна и относится к распределению значений), которая является максимальным числом сгенерированных вызовов пользователем и параметром масштаба. В частности, в моем случае это значение равно 500, что означает, что один пользователь может генерировать максимум 500 вызовов.

Маурицио

(часть 2/2) Другим параметром, который нужно установить, является параметр масштаба для Zipf, который определяет разброс распределения (в моем случае это значение равно 0,9). У меня есть все параметры (размер выборки, частота на пользователя и т. Д.), Но я должен убедиться, что мой набор данных соответствует распределению zipf.

Маурицио

так что вы, очевидно, перенормируете распределение на , поскольку для того, что я бы назвал «усеченным Zipf», параметр масштабирования 0,9 был бы невозможен , Если вы можете сгенерировать много этих данных и у вас «только» есть 500 возможных результатов, почему бы просто не использовать критерий соответствия хи-квадрат? Поскольку у вашего дистрибутива длинный хвост, вам может потребоваться довольно большой размер выборки. Но это был бы один из способов. Другим быстрым и грязным методом было бы убедиться, что вы получаете правильное эмпирическое распределение для небольших значений числа вызовов.

\sum_{i = 1}^{500} i^{- 0.9}

$\sum_{i=1}^{500} i^{-0.9}$

кардинал

Ответы:

Обновление : 7 апреля 2011 г. Этот ответ становится довольно длинным и охватывает несколько аспектов проблемы. Однако до сих пор я сопротивлялся, разбивая его на отдельные ответы.

В самом низу я добавил обсуждение производительности Пирсона для этого примера. $\chi^2$

Брюс М. Хилл написал, пожалуй, «основополагающую» статью об оценке в Zipf-подобном контексте. В середине 1970-х он написал несколько статей на эту тему. Тем не менее, «оценщик Хилла» (как он теперь называется) по существу полагается на статистику максимального порядка выборки, и поэтому, в зависимости от типа присутствующего усечения, это может привести к некоторым проблемам.

Основная статья:

Б. М. Хилл, простой общий подход к выводу о хвосте распределения , Ann. Стат. 1975 г.

Если ваши данные изначально являются Zipf, а затем усечены, то хорошее соответствие между распределением степеней и графиком Zipf может быть использовано в ваших интересах.

В частности, распределение степеней - это просто эмпирическое распределение числа раз, когда просматривается каждый целочисленный ответ,

d_{i} = \frac{# {j : X_{j} = i}}{n} .

$d_i = \frac{\#\{j: X_j = i\}}{n} .$

Если мы построим это против на графике log-log, мы получим линейный тренд с наклоном, соответствующим коэффициенту масштабирования. $i$

С другой стороны, если мы строим график Zipf , где мы сортируем выборку от наибольшего к наименьшему, а затем наносим значения на их ранги, мы получаем другой линейный тренд с другим наклоном. Однако склоны связаны между собой.

Если является коэффициентом закона масштабирования для распределения Zipf, то наклон на первом графике равен а наклон на втором графике равен . Ниже приведен пример графика для и . Левая панель - это распределение степеней, а наклон красной линии равен . Правая часть - это график Зипфа, с наложенной красной линией, имеющей наклон . $\alpha$ $-\alpha$ $-1/(\alpha-1)$ $\alpha = 2$ $n = 10^6$ $-2$ $-1/(2-1) = -1$

Графики распределения степени (слева) и Zipf (справа) для выборки iid из распределения Zipf.

Таким образом, если ваши данные были усечены, так что вы не видите значений, превышающих некоторый порог , но в противном случае данные распределяются по Zipf, а достаточно велик, тогда вы можете оценить по распределению степеней . Очень простой подход состоит в том, чтобы подогнать линию к графику log-log и использовать соответствующий коэффициент. $\tau$ $\tau$ $\alpha$

Если ваши данные усечены так, что вы не видите маленьких значений (например, как много фильтрации выполняется для больших наборов веб-данных), то вы можете использовать график Zipf для оценки наклона в масштабе log-log и затем " отступить "показатель масштабирования. Скажите, что ваша оценка наклона на графике Zipf равна . Тогда одной простой оценкой коэффициента закона масштабирования является $\hat{\beta}$

\hat{α} = 1 - \frac{1}{\hat{β}} .

$\hat{\alpha} = 1 - \frac{1}{\hat{\beta}} .$

@csgillespie дал одну недавнюю статью, написанную в соавторстве с Марком Ньюманом в Мичигане на эту тему. Похоже, он публикует много похожих статей на эту тему. Ниже приведено еще несколько ссылок, которые могут представлять интерес. Ньюман иногда не делает ничего разумного статистически, поэтому будьте осторожны.

MEJ Newman, Степенные законы, распределения Парето и закон Ципфа , Современная физика 46, 2005, с. 323-351.

Mitzenmacher, краткая история генеративных моделей для степенного закона и логнормальных распределений , Internet Math. том 1, нет. 2, 2003, с. 226-251.

К. Найт, простая модификация оценки Хилла с приложениями к устойчивости и снижению смещения , 2010.

Приложение :

Вот простая симуляция в чтобы продемонстрировать, что вы можете ожидать, если вы взяли выборку размером из своего дистрибутива (как описано в вашем комментарии ниже вашего исходного вопроса). $R$ $10^5$

> x <- (1:500)^(-0.9)
> p <- x / sum(x)
> y <- sample(length(p), size=100000, repl=TRUE, prob=p)
> tab <- table(y)
> plot( 1:500, tab/sum(tab), log="xy", pch=20, 
        main="'Truncated' Zipf simulation (truncated at i=500)",
        xlab="Response", ylab="Probability" )
> lines(p, col="red", lwd=2)

Получившийся сюжет

«Усеченный» график Zipf (усеченный при i = 500)

Из графика видно, что относительная погрешность распределения степеней для (или около того) очень хорошая. Вы могли бы сделать формальный критерий хи-квадрат, но это не строго сказать вам , что данные следуют предуказанных распределения. Это только говорит вам, что у вас нет доказательств, чтобы сделать вывод, что они этого не делают . $i \leq 30$

Тем не менее, с практической точки зрения, такой сюжет должен быть относительно убедительным.

Приложение 2 : Давайте рассмотрим пример, который Маурицио использует в своих комментариях ниже. Предположим, что и , с усеченным распределением Zipf, имеющим максимальное значение . $\alpha = 2$ $n = 300\,000$ $x_{\mathrm{max}} = 500$

Мы рассчитаем статистику Пирсона двумя способами. Стандартный способ - через статистику где - это наблюдаемые значения значения в образце и . $\chi^2$

X^{2} = \sum_{i = 1}^{500} \frac{(O_{i} - E_{i})^{2}}{E_{i}}

$X^2 = \sum_{i=1}^{500} \frac{(O_i - E_i)^2}{E_i}$

O_{i}

$O_i$

i

$i$

E_{i} = n p_{i} = n i^{- α} / \sum_{j = 1}^{500} j^{- α}

$E_i = n p_i = n i^{-\alpha} / \sum_{j=1}^{500} j^{-\alpha}$

Мы также вычислим вторую статистику, сформированную первым объединением счетчиков в ячейках размером 40, как показано в электронной таблице Маурицио (последняя ячейка содержит только сумму из двадцати отдельных конечных значений.

Давайте нарисуем 5000 отдельных выборок размера из этого распределения и вычислим используя эти две разные статистики. $n$ $p$

Гистограммы значений приведены ниже и выглядят достаточно однородными. Эмпирические коэффициенты ошибок типа I составляют соответственно 0,0716 (стандартный метод без объединения) и 0,0502 (метод с сортировкой), и ни один из них статистически значимо не отличается от целевого значения 0,05 для размера выборки 5000, который мы выбрали. $p$

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

Вот код $R$

# Chi-square testing of the truncated Zipf.

a <- 2
n <- 300000
xmax <- 500

nreps <- 5000

zipf.chisq.test <- function(n, a=0.9, xmax=500, bin.size = 40)
{
  # Make the probability vector
  x <- (1:xmax)^(-a)
  p <- x / sum(x)

  # Do the sampling
  y <- sample(length(p), size=n, repl=TRUE, prob=p)

  # Use tabulate, NOT table!
  tab <- tabulate(y,xmax)

  # unbinned chi-square stat and p-value
  discrepancy <- (tab-n*p)^2/(n*p)
  chi.stat <- sum(discrepancy)
  p.val    <- pchisq(chi.stat, df=xmax-1, lower.tail = FALSE)

  # binned chi-square stat and p-value
  bins <- seq(bin.size,xmax,by=bin.size)
  if( bins[length(bins)] != xmax )
    bins <- c(bins, xmax)

  tab.bin  <- cumsum(tab)[bins]
  tab.bin <- c(tab.bin[1], diff(tab.bin))

  prob.bin <- cumsum(p)[bins] 
  prob.bin <- c(prob.bin[1], diff(prob.bin))

  disc.bin <- (tab.bin - n*prob.bin)^2/(n * prob.bin)
  chi.stat.bin <- sum(disc.bin)
  p.val.bin <- pchisq(chi.stat.bin, df=length(tab.bin)-1, lower.tail = FALSE)

  # Return the binned and unbineed p-values
  c(p.val, p.val.bin, chi.stat, chi.stat.bin)
}

set.seed( .Random.seed[2] )

all <- replicate(nreps, zipf.chisq.test(n, a, xmax))

par(mfrow=c(2,1))
hist( all[1,], breaks=20, col="darkgrey", border="white",
      main="Histogram of unbinned chi-square p-values", xlab="p-value")
hist( all[2,], breaks=20, col="darkgrey", border="white",
      main="Histogram of binned chi-square p-values", xlab="p-value" )

type.one.error <- rowMeans( all[1:2,] < 0.05 )

кардинальный
источник

+1, отличный ответ как обычно. Вы должны назначить себя модератором, осталось еще 1 час :)

mpiktas

@mpiktas, спасибо за комплименты и поддержку. Я не уверен, что смог бы оправдать выдвижение своей кандидатуры уже и без того сильным списком кандидатов, которые, как правило, участвовали более активно и дольше, чем я.

кардинал

@cardinal, вот некоторые ссылки на Альтернативу оценщика Хилла: оригинал статья по Paulauskas и прослеживание по Vaiciulis и Gadeikis и Paulauskas . Эта оценка предположительно имела лучшие свойства, чем первоначальная оценка Хилла.

mpiktas

@mpiktas, спасибо за ссылки. Существует довольно много «новых и улучшенных» версий оценщика Хилла. Основным недостатком оригинального подхода является то, что он требует выбора «отсечки», где следует прекратить усреднение. Я думаю, что в основном это было сделано «наглазником», который открывает обвинения в субъективности. Одна из книг Резника о длиннохвостых дистрибутивах обсуждает это в некоторых деталях, если я помню. Я думаю, что это его более свежий.

кардинал

@cardinal, большое спасибо, вы очень добры и очень подробны! Ваш пример в R был очень полезен для меня, но как я могу выполнить формальный тест хи-квадрат в этом случае? (Я использовал тест хи-квадрат с другими распределениями, такими как равномерное, экспоненциальное, нормальное, но у меня много сомнений по поводу zipf .. Извините, но это мой первый подход к этим темам). Вопрос к модераторам: должен ли я написать еще один вопрос и ответ, например "Как выполнить тест хи-квадрат для усеченного распределения zipf?" или продолжить в этом Q & A возможно обновление тегов и заголовка?

Маурицио

Бумага

Clauset, A и др , степенные распределения в эмпирических данных . 2009

содержит очень хорошее описание того, как подходить к модели степенного закона. На соответствующей веб-странице есть примеры кода. К сожалению, он не дает код для усеченных дистрибутивов, но может дать вам указатель.

Кроме того, в статье обсуждается тот факт, что многие "степенные наборы данных" могут быть смоделированы одинаково хорошо (а в некоторых случаях лучше) с помощью нормального или экспоненциального распределения Log!

csgillespie
источник

К сожалению, в этой статье ничего не говорится о усеченном дистрибутиве. Я нашел несколько пакетов в R, которые просто имеют дело с параметром оценки Zipf (zipfR, VGAM), но усеченному дистрибутиву нужна «специальная обработка». В последнем предложении вы имели в виду, что можно моделировать степенной набор данных с, например, экспоненциальным распределением, а затем применить некоторый процесс оценки параметров для «усеченного» экспоненциального распределения? Я очень новичок в этой теме!

Маурицио

В статье авторы повторно анализируют различные наборы данных, где был установлен степенной закон. Авторы указывают, что в ряде случаев степенная модель не так уж велика, и альтернативное распределение было бы лучше.

csgillespie

После подробного ответа пользователя кардинала я выполнил тест хи-квадрат на моем предположительном усеченном распределении zipf. Результаты теста хи-квадрат представлены в следующей таблице:

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

Где StartInterval и EndInterval представляют, например, диапазон вызовов, а Observed - это количество вызывающих абонентов, генерирующих от 0 до 19 вызовов и т. Д. Тест хи-квадрат хорош до достижения последних столбцов, они увеличивают окончательный расчет, иначе до этого момента «частичное» значение хи-квадрат было приемлемым!

С другими тестами результат тот же, последний столбец (или последние 2 столбца) всегда увеличивает конечное значение, и я не знаю почему, и я не знаю, используют ли (и как) другой проверочный тест.

PS: для полноты, чтобы рассчитать ожидаемые значения ( Ожидаемые ), я следую совету кардинала следующим образом:

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

где X_i «s используются для расчета: x <- (1:n)^-S, то P_i » S , чтобы вычислить p <- x / sum(x)и , наконец, e_i (ожидаемый Н.Р. пользователей для каждого пг вызовов) получаетсяP_i * Total_Caller_Observed

а со степенью свободы = 13 добродетель хи-квадрата всегда отвергает гипотезу о том, что набор выборок следует распределению Zipf, потому что статистика тестов (в данном случае 64,14) больше, чем в таблицах хи-квадрат, "demerit" для последнего столбца. Графический результат сообщается здесь: введите описание изображения здесь

хотя точка усечения установлена на 500, максимальное значение получается равным 294. Я думаю, что окончательная «дисперсия» является причиной провала теста хи-квадрат.

ОБНОВИТЬ!!

Я пытаюсь выполнить тест хи-квадрат на предполагаемой выборке данных zipf, сгенерированной с помощью кода R, о котором сообщалось в ответе выше.

> x <- (1:500)^(-2)
> p <- x / sum(x)
> y <- sample(length(p), size=300000, repl=TRUE, prob=p)
> tab <- table(y)
> length(tab)
[1] 438
> plot( 1:438, tab/sum(tab), log="xy", pch=20, 
        main="'Truncated' Zipf simulation (truncated at i=500)",
        xlab="Response", ylab="Probability" )
> lines(p, col="red", lwd=2)

Связанный сюжет следующий: введите описание изображения здесь

Результаты теста хи-квадрат представлены на следующем рисунке: введите описание изображения здесь

и статистика теста хи-квадрат (44,57) слишком высока для проверки с выбранной степенью свободы. Также в этом случае окончательная «дисперсия» данных является причиной высокого значения хи-квадрат. Но есть процедура для проверки этого распределения zipf (независимо от моего «неправильного» генератора, я хочу сосредоточиться на образце данных R) ???

Maurizio
источник

@Maurizio, по некоторым причинам, я пропустил этот пост до сих пор. Можно ли как-нибудь отредактировать его и добавить сюжет, похожий на последний в моем посте, но используя ваши данные наблюдений? Это может помочь диагностировать проблему. Я думаю, что я видел еще один ваш вопрос, где у вас были проблемы с созданием равномерного распределения, так что, возможно, это относится и к этим анализам. (?) С уважением.

кардинал

@cardinal, я обновил результаты! Что вы думаете? Вопрос о равномерном распределении - это еще одна вещь, которую я должен уточнить, и я сделаю это сегодня или завтра;)

Маурицио

@ Маурицио, они были созданы случайно? Был ли ваш масштабный параметр как и раньше? Я использовал размер выборки 8454 и точку усечения 500 и сгенерировал 10000 таких выборок. Из этих 10000 максимальное наблюдаемое значение в выборке составило 500 для 9658 испытаний, 499 для 324 испытаний, 498 для 16 испытаний и 497 для 2 испытаний. Исходя из этого, я думаю, что что-то не так с вашей процедурой генерации. Если вы не использовали другой параметр масштаба.

S = 0.9

$S = 0.9$

кардинал

@ Маурицио, чтобы объяснить результаты, которые я выложил, учтите, что . Так, при размере выборки ожидаемое количество результатов со значением 500 составляет . Вероятность увидеть хотя бы один такой исход составляет . Обратите внимание, насколько близко это соответствует моделированию выше.

p = P (X_{i} = 500) \approx 4.05 \times 10^{- 4}

$p = \mathbb{P}(X_i = 500) \approx 4.05 \times 10^{-4}$

n = 8454

$n = 8454$

8454 \cdot 4.05 \cdot 10^{- 4} \approx 3.43

$8454 \cdot 4.05 \cdot 10^{-4} \approx 3.43$

1 - (1 - 0.000405)^{8454} \approx 0.9675

$1 - (1 - 0.000405)^{8454} \approx 0.9675$

кардинал

@cardinal, я также думаю, что в процедуре генерации есть что-то «неправильное» (моя цель - проверить, действительно ли этот генератор следует распределению Zipf). Я должен поговорить с дизайнерами проекта в эти дни.

Маурицио