СТАТИСТИЧЕСКИЕ СВОЙСТВА ПОКАЗАТЕЛЯ РАСПРЕДЕЛЕНИЯ ПАРЕТО Текст научной статьи по специальности «Математика»

Cloud of Science. 2020. T. 7. № 3 http:/ / cloudofscience.ru

Статистические свойства показателя распределения Парето

И. С. Пулькин, А. В. Татаринцев

МИРЭА - Российский технологический университет, 119454, Москва, пр-т Вернадского, 78

e-mail: pulkin@mirea.ru, tatarintsev@mirea.ru

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

Ключевые слова: показатель распределения Парето, функция распределения, моменты, числовые характеристики распределения.

1. Введение

Распределение Парето задается плотностью вероятностей

при х >0. При х <0 плотность вероятности и функция распределения равны нулю.

Это распределение впервые появилось в работах Вильфредо Парето по экономической статистике. Оно описывает распределение доходов в обществе, при условии, что все эти доходы превышают некоторый базовый уровень 8.

Впоследствии были открыты многие новые важные области применения этого распределения: лингвистика (частота употребления слов); сейсмология (закон Гу-

и функцией распределения

тенберга — Рихтера); распределение длин файлов в интернет-трафике; популярность имен; размер паводка; природные и техногенные катастрофы.

Степенные законы и связанные с ними распределения широко исследуются в последнее время. Многие области исследований и применений описаны, например, в [1].

Возможно, самой практически важной проблемой является исследование катастроф. Предсказание интенсивности катастроф, как природных, так и техногенных, — это задача, важность которой переоценить невозможно. Вот, например, что сказано по этому поводу в [2]: «...В ряду ущербов от катастроф изредка встречаются суперэкстремальные значения, несоизмеримые по величине со значениями для подавляющей части событий. Ущерб от этих суперэкстремальных событий сравним с суммарным ущербом от всех катастроф за тот же период времени».

2. Обсуждение проблемы и предварительные исследования

При планировании необходимо знать, какой интенсивности катастроф следует опасаться, а какие маловероятны. Для этого надо уметь определять параметры распределения Парето по предыдущим наблюдениям. Это приводит к такой математической постановке задачи: получена выборка х1 ,...хп, подчиняющаяся распределению Парето. Требуется определить или с достаточной точностью оценить параметры а и 9.

Сразу заметим, что достаточно определить параметр а, поскольку параметр 9 обычно известен заранее: это — минимальный порог, и значения, меньшие 9, просто не попадают в выборку. Если исследуем, скажем, землетрясения, то в выборку попадают только те из них, чья интенсивность (энергия, магнитуда) больше определенного порога.

Если возьмем другую выборку, с другим порогом, то параметр а окажется тем же самым. Это следует из такого свойства распределения Парето: если X — случайная величина, распределенная по Парето с параметрами а и 9, то условное распределение X при условии X > х также подчиняется распределению Парето с параметрами а и х •

Таким образом, следует сосредоточиться на задаче оценивания по выборке параметра а, считая, что параметр 9 известен.

Для оценки параметров распределения по выборке наиболее часто используется метод максимального правдоподобия. В случае распределения Парето плотность вероятности того, что элементы выборки примут значения х,..., хп, если параметры распределения Парето равны а и 9, равна

0(х1,...хп |а,0) = | а

п ( 0п \

Vх ' .••' хп у

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

п

а = -

1пх + ... + 1пхя -п 1п0 Здесь х1,..., хп — выборка; п — объем выборки.

Эту формулу можно представить в виде

п

а =-,

1п( х1/ 0) + ... + 1п( хя/ 0)

откуда следует, что без ограничения общности можно считать, что 0 = 1.

Настоящая работа посвящена исследованию статистических свойств оценки максимального правдоподобия параметра а. В [2] обсуждаются различные, в основном эмпирические, подходы к этой оценке, в частности, предпринимается попытка интервального оценивания. Другие подходы см. например в [3] и [4].

В [5] описан следующий вычислительный эксперимент. Генерировалось 5000 выборок, подчиненных распределению Парето, с параметрами а = 0.8 и а = 2.5. Объемы выборок были равны 10, 20 и 50. Для каждой выборки вычислялась оценка а. Результаты вычислительного эксперимента указывают на то, что оценка а — смещенная и позволяют предложить эмпирическую формулу для несмещенной оценки:

п -1 ,

а =-а.

п

Результаты численного эксперимента можно представить в табл. 1.

Таблица 1. Результаты численного эксперимента оценки максимального правдоподобия

а = 0.8

Объем выборки 10 20 50

ОМП а 0.890 0.843 0.813

Исправленная ОМП а 0.801 0.801 0.797

а = 2.5

Объем выборки 10 20 50

ОМП а 2.782 2.635 2.550

Исправленная ОМП а 2.504 2.503 2.499

В табл 1. указана как оценка максимального правдоподобия (ОМП) а, так и исправленнная в соответствии с предыдущей формулой оценка а.

Если брать различные выборки х, . ••, X, то несмещенная оценка максимального правдоподобия становится случайной величиной. Аналитически свойства этой случайной величины не исследовались, в частности, никакие ее свойства не нашли отражения в подробном и всеобъемлющем справочнике [6]. Исследованию свойств этой случайной величины и посвящена настоящая статья.

Некоторые свойства этого распределения были установлены авторами в [7]. В частности, там строго аналитически доказано, что оценка а — смещенная, и вместо нее следует рассматривать оценку

а = а(я-1)/я.

В [7] доказано, что это действительно несмещенная оценка, то есть ее математическое ожидание совпадает с истинным значением параметра а. При этом обе оценки состоятельные. Кроме того, в [7] авторами была доказана следующая теорема.

Пусть х1,..хп — выборка, подчиняющаяся распределению Парето с параметрами а, 9 и объемом выборки п. Случайная величина, равная исправленной оценке максимального правдоподобия для параметра распределения а

х =-П—1-

1пх + ••• + 1пх —п 1п9

имеет функцию распределения

п—1 1

^(х) = Р(X < х) = е—а7х (а/х)к

к=0 к!

и плотность вероятности

/ Лп+1 1 ( а | ~—а / х

) аГ(п) V х )

В этих формулах и введен новый параметр а = а9(п — 1), объединяющий параметры исходного распределения Парето.

Далее устанавливаются основные свойства этого распределения.

На рис. 1 и 2 представлены графики плотности вероятности этого распределения для а = 2.5 и а = 0.8 при различных значениях п, при этом 9 = 1. Как ранее указывалось, всегда можно считать, что 0 = 1, поскольку общий случай приводится к этому случаю заменой х; =х;/9. В статье, однако с целью сохранения общности, рассматривается случай произвольного 9.

Математическое ожидание для этого распределения равно МХ = а 9, дисперсия

ох=-(а9)!

п — 2

На рис. 3 и 4 представлены функции распределения при различных значениях а и п.

Рисунок 1. Графики плотности вероятностей распределения с а = 2.5 (красный) и а = 0.8 (зеленый), п = 20

Рисунок 2. Графики плотности вероятностей распределения с а = 2.5 (красный)

Рисунок 3. Графики функции распределения Рисунок 4. Графики функции распределения с а = 2.5 (красный) и а = °.8 (зеленый), с а = 2.5 (красный) и а = 0.8 (зеленый),

п = 20 п = 50

3. Основные результаты

Вычислим значения числовых характеристик распределения, полученного в [7].

Мода распределения, точка в которой плотность вероятности максимальна, получается из условия р '(х) = 0 и равна:

а п —1 „

хо =-7 =-7'а9-

п +1 п +1

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

Для вычисления медианы распределения заметим, что функция распределения зависит только от двух параметров, а/х и п. Обозначим у = а/х и найдем медиан-

ное значение ут, при котором функция распределения равна 1/2. После этого медиана находится по формуле

Хт =«% " 1)/Ут-

Эта медиана, разумеется, зависит от п. График этой зависимости, полученный в результате численных расчетов, приведен на рис. 5. При малых значениях п график малоинформативен, а с ростом п медиана быстро сходится к значению а.

Рисунок 5. Зависимость моды (синим) и медианы (фиолетовым) распределения от п, а = 2.5

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

Рассмотрим асимптотику плотности вероятности при больших п. Для этого выполним предельный переход в данном выражении следующим образом. Плотность вероятности р(х) представим в виде экспоненциальной функции:

1 „ ^

р( х) = —ехр (-ф( х) ); ф( х) = - + (п + 1)1п аГ(п) х

х V а

Легко проверить, что функция ф(х) имеет точку минимума в моде х0. При увеличении расстояния от моды до точки х плотность вероятности быстро спадает, так как ф(х) при (х - х0) или х ^+0. Скорость стремления к беско-

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

ф(х) = ф(х0) + ф'(х0)(х - х0) +

ф''( хо) 2!

(х - х0)2+ О ((х - х0)2).

Так как ф'( х0) = 0, то второе слагаемое в разложении отсутствует. Ограничившись квадратичным по (х - х0) слагаемым в разложении функции ф(х) мы получим приближенное, с точностью до экспоненциальных множителей ехр(-Сп); С > 0, асимптотическое распределение имеющее гауссову форму:

1 -'"Ч) _____( (х хо) ^

Рат (х) = ^7Т' е хо)' ехр

аГ (п)

2а2

Квадратичное отклонение этого распределения зависит от объема выборки и

других параметров исходного распределения Парето:

а=^ а9)

и стремится к нулю при п ^ +да. Предэкспоненциальный нормировочный множитель получившегося распределения при больших п эквиалентен нормировочному коэффициенту нормального распределения:

1 -ФСх ) 1 (п + 1Т+1 1

N =--е ф(хо) =--1 - I ~ ,_ при п^+да.

аГ(п) аГ(п) ^ е ) у2л а

Дисперсия асимптотического распределения также эквивалентна дисперсии исходного распределения при больших п, хотя имеются и некоторые различия в следующих членах разложения по 1/ п:

Оа8утр(Х) = а2= ' (а9)2~ Б(Х) = ^ п -++да

(п +1)3 и - 2

Таким образом, получим, что с точностью до экспоненциальных бесконечно малых величин при больших п имеет место распределение с плотностью вероятности:

Ра8утр (х) = /Г1 еХР (— (х — х0 )7(2а2 )). а

Это распределение удовлетворяет условию нормировки на всей числовой оси х е (—да; +да). Исходное же распределение вероятности имеет плотность отличную от нуля на интервале х е [0; +да). Такое расширение области определения допустимо, так как изменяет плотность вероятности на экспоненциальную бесконечно малую величину. Так как нормальное распределение симметрично, то мода и медиана совпадают с точностью до экспоненциальных бесконечно малых при больших п.

Математическое ожидание случайной величины X и ее квадрата, как уже было установлено в [7], равны:

а

М( X) =-= а9;

п — 1

М( X2) =---= —' (а9)2.

(п — 1)(п — 2) п — 2

Продолжая вычисление моментов случайной величины, легко получить,

а3 (п-1)2

М ( X 3) =-^-= --(а9)3;

' (п — 1)(п — 2)(п — 3) (п — 2)(п — 3)

М(Х4) =---=---(аЭ)4.

(п -1)(п - 2)(п - 3)(п - 4) (п - 2)(п - 3)(п - 4)

В общем случае, момент произвольного порядка для случайной величины X также легко вычисляется:

М(X-) = .ар = (п-1)Р-Г(п-р) .(ае)р.

4 ' Г (п) Г (п)

Для старших центральных моментов распределения имеем выражения:

4 • (аЭ)3

= М (X - М (X ))3 =

(п - 2)(п - 3):

* = м (х - М(х))4 =. (3п+15) •(ае)4

(п - 2)(п - 3)(п - 4) Коэффициент асимметрии:

4^ п - 2

73 = Б(X)3/2 " (п - 3) .

Коэффициент эксцесса:

ц4 30п - 66

7 =-4--3 =-

Б(Х )2 (п - 3)(п - 4)

положительный, так как п > 4.

Дифференциальная энтропия определяется с точностью до постоянного слагаемого, поэтому можно отбросить постоянную часть нормировочного множителя под знаком логарифма, поскольку этот сомножитель даст слагаемое, не зависящее от п. Для нашего распределения получим:

1 +да ( \п+1

Н(X) =-| I - I /х (1п ((п - 1)Г(п)) - (п +1) 1п (а/х) + (-х)) с1х .

аГ(п) 0 ^ х)

Снова сделаем замену переменной интегрирования а/х = у, получим:

л +да

Н(X) = од 1 У" 1е-у • (1п ((п - 1)Г(п)) + у - (п + 1)1п у)йу =

-I +да

= 1п ((п - 1)Г(п)) + п - п— | уп-1 е-у • 1п у йу = Г (п) 0

= 1п ( (п - 1)Г(п)) + п - п+1 • Г '(п).

Используя свойства гамма-функции Эйлера, описанные в [8], можно записать выражение для дифференциальной энтропии в виде:

Н(X) = 1п ((п - 1)Г(п)) + п - (п + 1)у(п)

или немного по другому

Н (X) = 1п (Г(п +1) — Г(п)) + п — (п +1)у(п). Здесь используются обозначения

, , Г '(х)

X) = -Сх) =—1п Г( х) Г(х) ах

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

п—1 2

у(п) = —у + ^ _ п > 2,

-1 к

где у = 0.5772... — постоянная Эйлера.

На рис. 6 приведен график дифференциальной энтропии в зависимости от п.

Рисунок 6. Графики дифференциальной энтропии в зависимости от объема выборки п

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

(4) = | егхр( х) —х.

Поскольку плотность вероятности имеет степенное убывание при х ^ +да, любой возрастающий экспоненциальный множитель приводит к расходящимуся интегралу. Поэтому производящая функция определена при 4 < 0. Моменты случайной величины выражаются как

—1 —хп

+да / \ п+1

М (Xn) = ^X (')1

^X ($) = I ' е*— а 7 х—х.

аГ(п) 0 V х )

0

После замены а/х = у интеграл приводится к виду:

1

Рх С) = ^ | У""1- е-у+а,/уОу. г(п) 0

Это интегральное представление для одной из модифицированных функций Бесселя:

Здесь Кп (•) — функция Макдональда (модифицированная функция Бесселя 3-го рода). По прежнему параметр равен а = а9(п -1) и г < 0.

Характеристическую функцию исследуемого распределения можно вычислить как мнимое продожение производящей функции (?) = (/?), то, тогда:

ф*с )=- Кп" ).

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

Кроме того, оказалось, что исследуемое распределение, подобно многим часто используемым распределениям, таким, например, как распределения Рэлея, Нака-муры и Вейбулла, тесно связано с гамма-распределением.

Вероятность попадания случайной величины X в интервал х < X < х2 задается интегралом:

X Х2 / \п+1

Р(X < х < х) = | р(X) сХ =-11 а I - е-а,ХсХ.

X аГ(п) х V х )

Заменив переменную у = а/х и пределы интегрирования у 2 = а/х 2, получим:

1 л л

Р(х <х < Х2) = — |у"-1 - е-уСу = | Рг(у)Су = Р(у2 < У < у).

( ) У2 У2

Здесь рг (у) — плотность вероятности гамма распределения

Рг(у) = ^у"-1 - е-".

Такая связь распределений говорит о родственном, но не тождественном отношении двух случайных величин. Таким образом, если случайная величина X имеет исследуемое нами распределение, то случайная величина У = а/Х имеет гамма распределение.

5. Заключение

Для решения задачи об определении показателя степенного закона по заданной выборке часто используется оценка максимального правдоподобия, которая в свою очередь приводит к случайной величине. Закон распределения этой величины которой был установлен в предыдущей публикации авторов [7].

В статье устанавливаются числовые характеристики этого распределения: мода, медиана, центральные моменты, а также дифференциальная энтропия и характеристическая функция. Получена аналитическа оценка моментов случайной величины как функций Макдональда, т. е. модифицированной функция Бесселя 3-го рода; а также характеристическая функция исследуемого распределения как мнимое продожение производящей функции моментов.

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

Литература

[1] Шредер М. Фракталы, хаос, степенные законы. Миниатюры из бесконечного рая. — Ижевск : НИЦ «Регулярная и хаотическая динамика», 2001.

[2] Владимиров В. А., Воробьев Ю. Л., Салов С. С. Управление риском: Риск. Устойчивое развитие. Синергетика / Препринт. — М. : ИПМ им. М. В. Келдыша, 2003. (http://www.keldysh.ru/papers/2003/source/book/gmalin/titul.htm)

[3] Белых В. В., Муравьев В. В., Степанов В. А. Применение обобщенного распределения Парето для определения структурных параметров графита в чугуне // Вестник ИжГТУ имени М. Т. Калашникова. 2019. Т. 22. № 1. С. 52-61.

[4] Волченкова И. В., Клебанов Л. Б. Характеризации распредения Парето свойствами соседних порядковых статистик // Записки научных семинаров ПОМИ. 2019. Т. 486. С. 63-70.

[5] Пулькин И. С. Методы оценки параметров распределения Парето // Математика, физика, информатика и их приложения в науке и образовании Международная школа-конференция молодых ученых: сборник тезисов докладов. — М. : МИРЭА, 2016. С. 123-124.

[6] Кобзарь А. И. Прикладная математическая статистика. Для инженеров и научных работников. — М.: ФИЗМАТЛИТ, 2005.

[7] Пулькин И. С., Татаринцев А. В. Свойства оценки максимального правдоподобия показателя распределения Парето // Российский технологический журнал. 2018. Т. 6. № 6. С. 74-83.

[8] Справочник по специальным функциям / Под ред. М.Абрамовица и И.Стигана. — М. : Наука, 1979.

Авторы:

Игорь Сергеевич Пулькин — кандидат физико-математических наук, доцент, доцент кафедры высшей математики, МИРЭА — Российский технологический университет

Андрей Владимирович Татаринцев — кандидат физико-математических наук, доцент, доцент кафедры высшей математики-2, МИРЭА — Российский технологический университет

Statistical properties of the Pareto distribution indicator

I. S. Pulkin, A. V. Tatarintsev

MIREA — Russian Technological University, 78, Vernadsky avenue Moscow, 119454 e-mail: pulkin@mirea.ru, tatarintsev@mirea.ru

Abstract. In this paper, we study the statistical properties of a random variable that is unbiased maximum likelihood estimation of the Pareto distribution indicator. Power distribution laws, such as the Pareto distribution, recently, researchers in various fields of science have been attracting close attention and technology, from Economics and linguistics to Internet traffic analysis. Therefore, the problem of determining the power law indicator for a given sample is exceptional practical importance. The most commonly used maximum likelihood estimation for such problems is results in a random value, the distribution law of which was established in the previous publication of the authors. This article sets the numerical characteristics of this distribution: mode, median, Central moments, as well as differential entropy and characteristic function. Keywords: index of the Pareto distribution, distribution function, moments, numerical characteristics of the distribution.

References

[1] Shreder M. (2001) Fraktaly, haos, stepennye zakony. Miniatyury iz beskonechnogo raya. [Rus]

[2] Vladimirov V. A., Vorob'ev YU. L., Salov S. S. (2003) Upravlenie riskom (Preprint, IPM) [Rus]

[3] Belyh V. V. et al. (2019) VestnikIzhGTUimeniM. T. Kalashnikova, 22(1):52-61. [Rus]

[4] Volchyonkova I. V., Klebanov L. B. (2019) Zapiski nauchnyh seminarov POMI, 486:63-70. [Rus]

[5] Pulkin I. S. (2016) Metody ocenki parametrov raspredeleniya Pareto. InMatematika, fizika, informatika i ih prilozheniya v nauke i obrazovanii Mezhdunarodnaya shkola-konferenciya molodyh uchenyh (MIREA), pp. 123-124. [Rus]

[6] Kobzar' A. I. (2005) Prikladnaya matematicheskaya statistika. Dlya inzhenerov i nauchnyh rabotnikov (FIZMATLIT). [Rus]

[7] Pulkin I. S., Tatarincev A. V. (2018) Rossijskij tekhnologicheskij zhurnal, 6(6):74-83. [Rus]

[8] AbramovicaM., Stigana I. (Eds.) (1979) Spravochnikpo special'nym funkciyam. (Nauka). [Rus]