Приветствую уважаемых читателей, случайно или нет наткнувшихся на эту статью. Адресована она всем тем, кто исследует распределение разнообразных эмпирических данных.
В моей публикации вас ждут: небольшой ликбез по теории вероятностей, развивающий её интуитивное и практическое понимание; детективная история о том, как решение дифференциального уравнения привело нас к двум важнейшим статистическим характеристикам выборочных последовательностей; ревизионизм в области проверки распределений на соответствие гауссовскому и объяснение пресловутого p-value с нуля за 5 минут (можете засекать).

Давайте начнем с ликбеза.

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

Сопутствующие определения и термины

Определение 1 Случайная величина X — переменная, принимающая, в зависимости от случая, те или иные числовые значения с определенными вероятностями.

Определение 2 Выборка (или статистический ряд) — последовательность из N наблюдений случайной величины X. Представляет собой набор [x_{1}, x_{2}, ... x_{n}], где — i-ый элемент выборки, а n — объем выборки.

Определение 3 Мода — самое часто встречающееся значение во множестве наблюдений случайной величины X. Название как бы намекает.

Определение 4 Математическое ожидание (выборочное среднее) — среднее арифметическое выборки, рассчитываемое по формуле:

\mu = \frac{1}{n} \sum_{i=1}^{n} x_i = \frac{1}{n}(x_{1} + ... + x_{n}).

Определение 5 Среднеквадратичное отклонение (стандартное отклонение) — показатель разброса значений случайной величины относительно её математического ожидания. Рассчитывается по формуле:

\sigma = \sqrt{\frac{1}{n}\sum_{i=1}^{n} (x_i - \mu) ^ 2}.

Определение 6 k-м центральным моментом случайной величины Xназывается величина

\begin{equation} \mu_{k} = \frac{1}{n} \sum_{i=1}^{n} (x_{i} - \mu) ^ {k}.\end{equation}

Нетрудно увидеть, что \mu_{2} есть не что иное, как квадрат среднеквадратичного отклонения \sigma ^ 2 и показывает разброс распределения вокруг среднего значения. А с моментами \mu_{3}и \mu_{4} мы еще пересечемся.

Определение 7 Распределение вероятностей — это закон, описывающий какие значения может принимать случайная величина и соответствующие вероятности появления этих значений.

Определение 8 Нормальное (гауссовское) распределение — это распределение вероятностей, которое играет важную роль во многих областях науки. Случайная величина подчиняется нормальному распределению, когда она подвержена влиянию огромного числа случайных помех. Такая ситуация крайне распространена, поэтому в природе это самое часто встречающееся распределение.

Нормальное распределение зависит от двух параметров — смещения (соответствует математическому ожиданию \mu) и разброса (соответствует стандартному отклонению \sigma).С математической точки зрения нормальное распределение является целым семейством распределений.

Определение 9 Стандартное нормальное распределение — нормальное распределение с математическим ожиданием \mu = 0 и стандартным отклонением \sigma = 1.

Кто эти ваши кривые Пирсона? 

Теперь начнем нашу небольшую детективную историю. Представьте, что вам на почту анонимно прислали письмо. Во вложениях у него архив, содержащий энное количество наборов из 256 вещественных чисел каждый. В письме вас просят найти, к какому типу вероятностого распределения относятся эти статистические ряды. К счастью, вам дали одну зацепку — кривые Пирсона. Google по этому запросу в первой строчке выдает статью [1]. Из неё мы узнаем, что кривые Пирсона - это семейство решений дифференциального уравнения

\frac{1}{y}\frac{dy}{dx} = \frac{-(x - a)}{b_{0} + b_{1}x + b_{2}x ^ 2}

где a, b_{0}, b_{1}, b_{2} — коэффициенты уравнения.

Более того, кривые Пирсона позволяют аппроксимировать практически любое известное распределение вероятностей. Иными словами, для каждого распределения найдется такая кривая Пирсона, которая по общим чертам будет походить на данное распределение.
Возникает вопрос: а каким же образом устанавливается связь между распределением и соответствующей ему кривой Пирсона?
Ответ: через коэффициенты a, b_{0}, b_{1}, b_{2} в уравнении выше, которые согласно [1] могут быть вычислены при помощи центральных моментов данного распределения из соотношений:

a = Mo, b_{0} = \frac{\mu_{2}(r + 1)}{r - 2}, b_{1} = -Mo = \frac{\mu_{3}(r + 2)}{2\mu_{2}(r - 2)}, b_{2} = -\frac{1}{r - 2},

где

r = \frac{6(\beta_{2} - \beta_{1} - 1)}{3\beta_{1} - 2\beta_{2} + 6}\beta_{1} = \frac{\mu_{3}^2}{\mu_{2} ^ 3}\beta_{2} = \frac{\mu_{4}}{\mu_{2} ^ 2}

\beta_{2} = \frac{\mu_{4}}{\mu_{2} ^ 2} — центральные моменты 2, 3 и 4 порядков соответственно, а Mo— мода.

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

Чтобы ответить на эти вопросы, полистаем статью [1] дальше. В ней вводится критерий\kappa(каппа Пирсона), который легко выражается через стандартизированные третий и четвертый моменты β_1 и β_2. Значения \kappa и вспомогательных параметров определяют 8 типов кривых Пирсона, соответствующих некоторым распространенным распределениям случайной величины X. В статье даже приводится график классификации областей в (β_1,β_2) координатах по этим 8 типам. Вы можете увидеть его на рис. 1:

Рис. 1. Графики для определения типа кривой Пирсона в зависимости от  и
Рис. 1. Графики для определения типа кривой Пирсона в зависимости от β_1 и β_2

Внимательный читатель заметил, что я писал о 8 типах, а на рис. 1 обозначено всего 7. Дело в том, что восьмой тип кривых Пирсона представляет собой точку с координатами (0, 3) и соответствует нормальному распределению. Иными словами, выборки, у которых β_1=0, а β_2=3 имеют нормальное распределение.

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

Рис. 2. Диаграмма рассеяния имеющихся статистических рядов. Красная прямая изображает верхнюю границу всех распределений. Синяя прямая — линию кривых III типа.
Рис. 2. Диаграмма рассеяния имеющихся статистических рядов. Красная прямая изображает верхнюю границу всех распределений. Синяя прямая — линию кривых III типа.

Нетрудно сравнить с рис. 1 и заметить, что все наши выборки находятся в области кривых I типа, причем в непосредственной близости к точке (0, 3). Ну и, естественно, теперь хочется понять, каким образом по высчитанным β_1 и β_2 можно оценить, насколько наше распределение соответствует нормальному.

Встречайте асимметрию и эксцесс

В широко известной в узких кругах статье "A Suggestion for Using Powerful and Informative Tests of Normality" [2] приводится подробное описание статистик \sqrt{β_1} (вместо β_1, с корнем формула симпатичнее выходит) и β_2, как мы знаем, задаваемые формулами (распишем чуть подробнее):

\sqrt{\beta_{1}} = \frac{\mu_{3}}{\mu_{2} ^ {3 / 2}} = \frac{\mu_{3}}{(\sigma ^ 2)^ {3 / 2}} = \frac{\mu_{3}}{\sigma ^ 3}

и

\beta_{2} = \frac{\mu_{4}}{\mu_{2} ^ 2} = \frac{\mu_{4}}{(\sigma ^ 2) ^ 2} = \frac{\mu_{4}}{\sigma ^ 4}

\beta_{2} = \frac{\mu_{4}}{\mu_{2} ^ 2} = \frac{\mu_{4}}{(\sigma ^ 2) ^ 2} = \frac{\mu_{4}}{\sigma ^ 4}де μ_3— 3-й центральный момент, μ_4 — 4-й центральный момент, а σ— среднее выборочное. Эти моменты, будучи стандартизированными (т.е. поделенными на σ в соответствующей степени), являются числовыми характеристиками асимметрии (skewness) распределения и его эксцесса (kurtosis), величины, показывающей, насколько ярко выражена вершина распределения в окрестности его среднего.

Преимущество этих статистик в их наглядности. Нормальное распределение симметрично, как мы знаем, поэтому для него \sqrt{\beta_{1}}=0. В то же время ненормальное асимметрично распределение имеет ненулевое значение асимметрии: \sqrt{\beta_{1}}>0 соответствует асимметрии вправо, а \sqrt{\beta_{1}}<0 соответствует асимметрии влево. Пусть вас не смущает, что арифметический корень может быть отрицательным, \sqrt{\beta_{1}}— это просто обозначение, если посмотрите на (9), то увидите, что в числите дроби третий момент который вполне может оказаться отрицательным числом. На рис. 3 приведена соответствующая иллюстрация.

Рис. 3. Иллюстрация распределений с разной асимметрией
Рис. 3. Иллюстрация распределений с разной асимметрией

Что касается четвертого момента, то у нормального распределения значение β_2 равно 3, а если быть точным 3(n-1)/(n+1), где n — объем выборки. Унимодальные распределения, т.е. распределения с одним пиком, чьи хвосты тяжелее (т.е. находятся ниже), чем у нормального, имеют β_2>3. Помимо этого, они обычно имеют более высокие пики в центре распределения. Унимодальные распределения с более легкими хвостами обычно имеют β_2<3. На рис. 4 проиллюстрированы описываемые распределения.

Рис. 4. Иллюстрация распределений с разной хвостатостью
Рис. 4. Иллюстрация распределений с разной хвостатостью

Теперь мы вполне понимаем смысл, стоящий за статистиками \sqrt{\beta_{1}} и β_2, и интуиция подсказывает нам, что немалую роль они играют в важном деле проверки распределений на нормальность. Но пока отложим этот вопрос и подумаем вот над чем.
Хорошо, рассчитали мы эти статистики для нашей выборки, получили числа. Вот, например, β_1=0.4. Вроде оно небольшое, близко к нули. А вроде и не такое уж и маленькое, вот 0.04 в 10 раз меньше. Понимаете проблему, да? На руках мы имеем абсолютные значения, а хотелось бы знать, если ли общепринятый порог, которые значения наших статистик могут перепрыгнуть ненароком. Вот тут нам и пригодится пресловутое p-value, ликбез которого был провокационно анонсирован уважаемым читателям в начале данной публикации. Ну что же, засекайте, товарищи.

Что такое p-value? 

Если вы посмотрите документацию по любому статистическому тесту на нормальность из библиотеки SciPy (например, тот же самый тест на асимметрию), то каждый из них возвращает ровно два числа: значение статистики (statistics) и загадочное число p-value. Ну с первым более-менее все ясно: почти каждый критерий нормальности определяется некоторой статистикой, и значение этой статистики по наблюдаемой выборке возвращается нам в качестве первого выходного значения функции.

А что же p-value? Тут нужно идти с самого начала. Тестирование выполняется путем сравнения значения статистики на конкретной выборки с нулевым распределением (в англоязычной литературе null distribution), то есть распределением уже самой статистики на множестве случайно сгенерированных выборок из нормального распределения. Другими словами, для того, чтобы построить нулевое распределение, нужно вычислить значения статистики на большом количестве статистических рядов из нормального распределения и построить из этих значений гистограмму. Важно отметить, что выборки для расчета нулевого распределения должны быть того же размера, что и наблюдаемый статистический ряд.

В уже упоминаемой выше статье [2], где представлен тест на асимметрию на основе нормированного третьего момента, вместо простого \sqrt{\beta_{1}} предлагается использовать в качестве статистики Z(\sqrt{\beta_{1}}) (в SciPy реализации используется этот вариант), где  Z-  довольно муторное отображение из R в R. Возникает 2 вопроса: что оно делает и зачем оно нужно? Давайте сравним нулевое распределение для \sqrt{\beta_{1}} и для Z(\sqrt{\beta_{1}}).
На рис. 5 продемонстрированы их гистограммы.

Рис. 5. Сравнение нулевых распределений для статистик  (custom) и  (scipy)
Рис. 5. Сравнение нулевых распределений для статистик \sqrt{\beta_{1}} (custom) и Z(\sqrt{\beta_{1}}) (scipy)

Гистограммы были получены с помощью теста Монте-Карло библиотеки SciPy, который случайным образом отбирает выборки из нормального распределения, рассчитывает нужную статистику на каждой из них и возвращает нулевое распределение. Несложно посчитать среднее значение и среднеквадратичное отклонение для каждого нулевого распределения. Для \sqrt{\beta_{1}} μ=0.001, σ=0.149, а для Z(\sqrt{\beta_{1}}) μ=0.002, σ=0.982. Видно, что нулевое распределение статистики из статьи [2] является стандартным нормальным распределением. А если почитать оригинальную статью [3], где впервые представлен skewness test, то можно понять, что это неслучайно: отображение Z есть не что иное, как общее преобразование Джонсона [4], которое приводит распределение подаваемой на вход статистики к стандартному нормальному.

Хорошо, с вопросом, что оно делает, мы разобрались. Остался второй вопрос: зачем оно нужно? Ответ прост: уже незачем. Вот у нас есть тест на нормальность на основе какой-то статистики. Мы можем вычислить её нулевое распределение, можем вычислить её значение применительно к наблюдаемой выборке. А как все-таки узнать по этому наблюдаемому значению статистики, нормальна наша выборка или нет? Для этого нам нужно оценить насколько наблюдаемое значение выбивается из общего распределения статистики для выборок из нормального распределения того же размера, иными словами, оценить долю значений нулевого распределения такой же или большей экстремальности, чем наблюдаемое значение. Эту долю и называют p-value.
Вообще говоря, существует 3 способа расчета p-value: подсчет доли значений статистики больше наблюдаемого (greater), меньше наблюдаемого (less) и больше наблюдаемого и меньше его значения, взятого со знаком минус (two-sided).
На рис. 6 приведена графическая интерпретация two-sided p-value.

Рис.6. Графическая интерпретация значения p-value
Рис.6. Графическая интерпретация значения p-value

Так вот, сейчас рассчитать эту долю не составляет труда, но в допотопные 60-80-ые, когда все эти методы тестирования на нормальность активно развивались, техника была ограниченная, поэтому рассчитать p-value для каждой степени уверенности, то есть доли более экстремальных значений, чем наблюдаемое, было трудоемко для всех известных распределений. Поэтому степени уверенности рассчитали довольно подробно только для стандартного нормального распределения, а статистики брали такими, чтобы их нулевое распределение соответствовало стандартному нормальному. Вот и весь секрет. Поэтому можно смело пользоваться статистиками \sqrt{\beta_{1}} и β_2 без дополнительных преобразований.

Практические результаты

Вообще говоря, для всех статистических тестов на нормальность можно улучшить точность вердикта. Дело в том, что при вычислении p-value никакого построения нулевого распределения не происходит — берется его приближение. Например, считается, как я уже говорил, что нулевое распределение у статистики Z(\sqrt{\beta_{1}}) похоже на стандартное нормальное, так же как и у Z(β_2). Или у не менее знаменитого критерия Д’Агостино Пирсона [5] cо статистикой Z(\sqrt{\beta_{1}})^2+Z(β_2 )^2 нулевое распределение близко к распределению χ^2 с двумя степенями свободы.
На рис. 7-8 приведены соответствующие гистограммы нулевого распределения вместе с их асимптотическими аппроксимациями.

Рис.7. Гистограмма нулевого распределения для
Рис.7. Гистограмма нулевого распределения для Z(β_2 )
Рис.8. Гистограмма нулевого распределения для
Рис.8. Гистограмма нулевого распределения для Z(\sqrt{\beta_{1}})^2+Z(β_2 )^2

Видно, что хоть аппроксимации эти и близки, но все же погрешность присутствует, особенно для выборок с малым количеством наблюдений. На рис. 9 приведен подтверждающий это пример. А значит, для того, чтобы повысить точность тестов, мы захотим использовать метод Монте Карло. И тут к нам на помощь и придет знание о том, что статистики \sqrt{\beta_{1}} и β_2 можно использовать без дополнительного преобразования Z, а это сакральное знание, в свою очередь, дает нам выигрыш в скорости. Давайте проверим, насколько он велик.

Рис. 9. Гистограмма нулевого распределения для статистики  с малым количеством экспериментов
Рис. 9. Гистограмма нулевого распределения для статистики Z(\sqrt{\beta_{1}})^2+Z(β_2 )^2 с малым количеством экспериментов

 Тест (статистика)

Время работы на одной выборке (в секундах)

Разница p-value на всех наборах (в %)

Our Skewness Test (\sqrt{\beta_{1}})

 1.9

 0.75

SciPy Skewness Test (Z(\sqrt{\beta_{1}})

3.0

Our Kurtosis Test (β_2)

 1.9

 0.16

SciPy Kurtosis Test (Z(β_2))

 3.1

Our Omnibus Test (β_1+β_2^2)

 3.7

 3.88

SciPy Normal Test (Z(\sqrt{\beta_{1}})^2+Z(β_2 )^2)

 6.4

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

Небольшое дополнение

В предыдущих двух главах ненароком упоминалось, что выявленные нами статистики лежат в основе критерия Д’Агостино Пирсона, настолько популярного, что в библиотеке SciPy его назвали normaltest. Но не только его. Наши \sqrt{\beta_{1}}и β_2 лежат, помимо прочего, в основе теста Харке-Бера, который использует оценку p-value для статистики

JB = n \Biggl( \frac{(\sqrt{\beta_{1}}) ^ 2}{6} + \frac{(\beta_{2} - 3) ^ 2}{24} \Biggr) (11)

Этот тест со звучным названием упоминается в статье [6] как самый мощный критерий для симметричных распределений с большим эксцессом (с тяжелыми хвостами). Ну и для фильтрации других типов ненормальных распределений он отлично подходит.

Итого, в данной статье мы 

  • познакомились с целой группой тестов на нормальность, называемой в англоязычной литература moment tests,

  • разобрались в том, что такое p-value и в общем принципе тестирования распределений на нормальность,

  • представили более быстрый способ применения вышеперечисленных критериев, чем тот, что реализован в SciPy.

References

[1]  Ralph B D’agostino, Albert Belanger, and Ralph B D’Agostino Jr.  A suggestion for using powerful and informative tests of normality.   The American Statistician, 44(4):316–321, 1990. 

[2]  Ralph B D’Agostino.  Transformation to normality of the null distribution of g1.   Biometrika, pages 679–681, 1970. 

[3]  Norman L Johnson.  Systems of frequency curves generated by methods of translation.   Biometrika, 36(1/2):149–176, 1949. 

[4]  RALPH D’agostino and Egon S Pearson.  Tests for departure from normality. empirical results for the distributions of b 2 and b 1.   Biometrika, 60(3):613–622, 1973. 

[5]  Bee Wah Yap and Chiaw Hock Sim.  Comparisons of various types of normality tests.   Journal of Statistical Computation and Simulation, 81(12):2141–2155, 2011.

Комментарии (2)


  1. kaptnemo
    25.06.2024 11:09

    "чьи хвосты тяжелее (т.е. находятся ниже)," -- "тяжёлые хвосты" распределения означают, что они выше, чем у нормального, а не ниже.


    1. zikmur Автор
      25.06.2024 11:09
      +1

      Дело, думаю, в определении. Не знаю, как в отечественной литературе, может названия прямо противоположные, я не проверял, если честно. Я пользовался терминологией из статьи [1]