Разберемся с измерением случайных величин (ч. 3/4)

Подбор подходящего распределения на основе экспериментальных данных

На английском такую операцию называют «distribution fitting«.

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

Статистические гипотезы

Статистической гипотезой называется любое предположение о виде или параметрах неизвестного закона распределения.

Проверяемую гипотезу обычно называют нулевой (или основной) и обозначают H_0. Наряду с нулевой гипотезой H_0 рассматривают альтернативную, или конкурирующую, гипотезу H_1, являющуюся логическим отрицанием H_0. Нулевая и альтернативная гипотезы представляют собой две возможности выбора, осуществляемого в задачах проверки статистических гипотез. Правило, по которому гипотеза H_0 отвергается или принимается, называется статистическим критерием или статистическим тестом.

Алгоритм выбора подходящего распределений

Существует множество способов (статистических критериев) определения подходящего распределения. Объединим их в общий алгоритм выбора распределений, начиная с самого простого и очевидного и заканчивая наиболее точным и надежным. В качестве экспериментальных данных будем использовать пример из первой части с распределением времен открытия веб-страниц.

Организационный момент: все вычисления будут проводится в веб-версии ipython – ipython notebook, которую вы можете запустить из консоли соответствующей командой: «ipython notebook». Это весьма удобная система для подобных задач, и, кстати, по общей концепции работы она похожа на рабочую тетрадь из дорогущей мат. системы Wolfram Mathematica, слава open-source!)

В конце я выложу файл с рабочей тетрадью (*.ipnb), в котором вы найдете все исходные коды из данной записи и сможете одним кликом запустить их у себя в своем ipython notebook.

1.1. Графический метод (ищем возможные подходящие распределения)

Суть графического метода в построении гистограммы полученных экспериментальных данных и последующим сравнении внешнего вида полученного графика с графиками плотности вероятностей (p.d.f.) различных законов распределения. И вправду, по сути гистограмма очень похожа на функцию плотности вероятности, только гистограмма весьма груба, дискретна и, значит, неточна, но для начала нам только это нужно для формирования гипотезы о действительном распределении нашей случайной величины.

Итак, запустим ipython с поддержкой построения графиков в самой рабочей тетради (что мега-круто!)) с помощью следующей команды в консоли:

ipython notebook --pylab inline

1. Получим необходимые данные с помощью следующего кода:

import httplib
import time

latencies = []
number_of_samples = 100
for i in range(number_of_samples):
    conn = httplib.HTTPConnection("livejournal.com")
    conn.request("GET", "/")
    start_time = time.time()
    r1 = conn.getresponse()
    latency = round((time.time() - start_time) * 1000)
    latencies.append(latency)
    print("{1}\\{2}. Задержка: {0}мс.".format(latency, i+1, number_of_samples))

2. Определим число участков разбиения и построим гистограмму времен открытия веб-страницы

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

k=\lfloor1+3.322*lg(n)\rfloor

либо просто

k=\sqrt{n},

где  k – число участков разбиения в гистограмме, n – число измерений.

В этот раз рассмотрим популярный сайт livejournal.com, он же ЖЖ.

import math
bins_number = round(10+math.sqrt(len(latencies)))
print(u"Число участков разбиения в гистограмме:{0}".format(bins_number))
<p>hist(latencies, bins_number)
title("Histogram of latencies of web-page load")
xlabel("Latency")
ylabel("Frequency")
show()

Получим следующее

3. Проанализируем гистограмму

Для начала отметим, что распределение несимметрично относительно значения с максимальной частотой появления (назовем такое значение «центром«). Правая часть распределения (или как говорят «хвост«) весьма растянута, т.е. вероятность большего значения случайного значения не так уж и мала. Такой хвост так и называют «длинный хвост» (long-tail). Если вероятность удаленного значения относительно велика, то такую ситуацию называют «тяжелый хвост» (heavy-tail).

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

4. Подберем параметры распределения и построим график распределения с наложением на гистограмму

И опять на помощь приходит numpy/scipy. В этой библиотеке есть порядка 80 распределений с функциями вероятности, плотности распределения и подбора параметров распределения на основе имеющихся данных.

Следующий код сделает ровно то, о чем гласит подзаголовок этого пункта.

import scipy.stats as ss
import scipy as sp
import numpy as np

hist(latencies, bins_number, normed=1)  # Построим нормализованную гистограмму

k = 1.5
print(ss.erlang.fit(latencies, k))
shape, loc, lambd = ss.erlang.fit(latencies, k)        # Оценим положение распределения

# Наберем больше промежуточных точек для плавности графика
x = np.linspace(np.min(latencies), np.max(latencies) + 10, number_of_samples * 2)

# Получим соответствующие значение по функции плотности вероятности для распределения Эрланга
distr_data = ss.erlang.pdf(x, k, loc=loc, scale=6)

# Построим график
plot(x, distr_data, 'r--')

Вы получите примерно следующее:

Левая хвост не особо вписался в наши данные, зато правый хвост подошел.

Любопытно, что если поменять параметр \lambda на 0.5, то мы получим распределение хи-квадрат с k-степенями свободы. Впрочем, это нам ничем здесь не поможет.)

В качестве альтернативного распределения выберем распределение Стьюдента, или t-распределение. По форме оно очень похоже на нормальное распределение и при значении единственного параметра, называемого степенью свободы (degree of freedom, df) близкого к 100-120 практически идентичного нормальному распределению. Главное его достоинство – более тяжелые хвосты, что полезно в ситуации малых выборок (меньше 30). На рисунке ниже голубая линия – нормальное распределение, красная – t-распределение.

Следующий код строит t-распределение для наших данных. Значения параметров распределения я подобрал вручную, поскольку функция fit() для t-распределения хорошо определяет только центр распределения, с другими параметрами у неё беда.

hist(latencies, bins_number, normed=1)  # Построим нормализованную гистограмму

df = 100

print(ss.t.fit(latencies, df))
tmp, loc, scale = ss.t.fit(latencies, df)        # Оценим положение распределения

# Наберем больше промежуточных точек для плавности графика
x = np.linspace(np.min(latencies)-10, np.max(latencies) + 10, number_of_samples * 2)

# Получим соответствующие значение по функции плотности вероятности для t-распределения
distr_data = ss.t.pdf(x, df, loc=loc, scale=8)

# Построим график
plot(x, distr_data, 'r--')

На глаз получилось не так уж плохо.)

Еще пачка важный обозначений

Перед рассказом о методе оценки коэффициента вариации введем пару обозначений, которые понадобятся нам в будущем.

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

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

Оценка – статистика, являющаяся основой для оценивания неизвестного параметра распределения.

Оценка может быть точечной и интервальной.

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

Интервальной называют оценку, которая определяется двумя числами – концами отрезка и вероятностью попадания истинного значения оцениваемого параметра в этот отрезок..

В таблице ниже приведены оценки различных параметров распределения.

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

Выборка (стат. совокупность)Случайная величина (генеральная совокупность)
Обозначение: {x_1,x_2,...,x_n}Обозначение: X
Характеристика: n – объем выборки, число элементов в выборкеХарактеристика: Закон распределения
Среднее арифметическоеArithmetic mean, mean, average

\overline{x}=\frac{1}{n}\sum_{i=1}^{n}x_i

Является оценкой мат. ожидания.

Математическое ожидание

Expected value, expectation, mathematical expectation, mean, first moment

\mu, E[X]

Cтандартное отклонение, стандартный разброс, выборочное стандартное отклонение
Sample standart deviation

S=\sqrt{\frac{1}{n-1}\sum_{i=1}^{n}(x_i-\overline{x})^2}

Является несмещённой оценкой стандартного отклонения случайной величины,т.е. мат. ожидание этой оценки равно оцениваемому параметру.

Стандартное отклонение

Standart deviation

\sigma=\sqrt{E[(X-\mu^2)]}

\sigma=\sqrt{E[X^2]-(E[X])^2}

Среднеквадратическое отклонение, среднеквадратичное отклонение, квадратичное отклонение
Standard deviation of the sample
S_x
S_N=\sqrt{\frac{1}{n}\sum_{i=1}^{n}(x_i-\overline{x})^2}

Также является оценкой стандартного отклонения распределения. В отличие от предыдущей оценки является смещенной оценкой.

Несмещенная (исправленная) выборочная дисперсия

Sample variance

S^2=\frac{1}{n-1}\sum_{i=1}^{n}(x_i-\overline{x})^2

Является несмещенной оценкой дисперсии.

Дисперсия

Variance

\sigma^2, \sigma_X^2, Var(X), D[X]

Выборочная дисперсия

Sample variance

S_n^2=\frac{1}{n}\sum_{i=1}^{n}(x_i-\overline{x})^2

Является смещенной оценкой дисперсии.

Среднее арифметического значение по полученным данным, среднее значение, выборочное среднее арифметическое, выборочное среднее – это простое среднее арифметическое, взятое из полученных данных. С терминологией здесь опять кто в лес, кто по дрова. Встречаются обозначения m*, \overline{M}, \overline{x}.

VarianceДисперсия – мера рассеивания значений относительно среднего значения. Малое значение дисперсии указывает на то, что разброс значений невелик относительно общего среднего значения.  Тут в литературе могут иметь в виду сразу два понятия:

Mean, expected value, expectationМатематическое ожидание – параметр закона распределения, по смыслу означает самое вероятное значение случайной величины. На графиках плотности распределения обычно на пике кривой. Есть у всех распределений (на сколько мне известно. Обозначается как M_x, m_x, M(x), \mu, E[X]. Считается для всех распределений по-разному. Для нормального распределения считается как среднее арифметическое по равномерно взятым значениям нормального распределения.

Фу-ух, с обозначениями пока всё. Теперь перейдем непосредственно к

1.2. Оценка коэффициента вариации

В случае, если вы уверены, что данные имеют, например, нормальное распределение, то в качестве простой альтернативы можно использовать метод оценки коэффициента вариации. Также его можно использовать дополнительно к графическому методу.

Коэффициент вариации V = \frac{\S_n}{\overline{x}} – это мера относительного разброса случайной величины; показывает, какую долю среднего значения этой величины составляет её средний разброс.

Полученное значение сравнивается со следующей таблицей и делается вывод, какое распределение лучше подходит:

ПРЕДЕЛЫ ИЗМЕНЕНИЯ КОЭФФИЦИЕНТА ВАРИАЦИИ VXЗАКОН РАСПРЕДЕЛЕНИЯ СЛУЧАЙНОЙ ВЕЛИЧИНЫ Х
Vx <= 0.3Нормальный
0.3 <= Vx < 0.4Гамма-распределение
0.4 <= Vx < 1Вейбулла
Vx = 1 Экспоненциальный

Попробуем метод на практике!

Тут я должен извиниться, я посеял данные с прошлого раза, поэтому собрал их снова, и теперь они выглядят немного по-другому:

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

Получилось даже лучше, чем раньше.)

Ну а теперь, наконец, посчитаем коэффициент вариации!

Данные получили следующие:

[189.0, 189.0,

186.0, 189.0, 189.0, 191.0, 190.0, 186.0, 188.0, 186.0, 187.0, 188.0,

186.0, 189.0, 208.0, 186.0, 194.0, 187.0, 186.0, 193.0, 190.0, 194.0,

186.0, 187.0, 201.0, 188.0, 188.0, 187.0, 186.0, 188.0, 189.0, 188.0,

185.0, 188.0, 191.0, 193.0, 194.0, 191.0, 188.0, 192.0, 187.0, 190.0,

189.0, 188.0, 193.0, 190.0, 188.0, 193.0, 187.0, 185.0, 186.0, 189.0,

190.0, 190.0, 189.0, 186.0, 186.0, 186.0, 187.0, 191.0, 186.0, 189.0,

188.0, 188.0, 194.0, 186.0, 195.0, 187.0, 185.0, 190.0, 187.0, 187.0,

188.0, 187.0, 188.0, 190.0, 188.0, 189.0, 189.0, 193.0, 187.0, 193.0,

186.0, 189.0, 188.0, 190.0, 186.0, 190.0, 188.0, 186.0, 185.0, 189.0,

189.0, 189.0, 194.0, 187.0, 192.0, 187.0, 190.0, 191.0]

Посчитаем среднее \overline{x}, стандартное отклонение S и коэф. вариации V_x

n = number_of_samples
mean = np.sum(latencies)/n
print("mean = {0}".format(mean))
S = np.sqrt(np.sum((latencies-np.mean(latencies))**2)/n-1)
print("S = {0}".format(S))
V_x = S/mean
print("V_x = {0}".format(V_x))

Получим \overline{x} = 189.01, S = 3.16384260038 и V_x = 0.016739022276

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

Критерии согласия

Критерий согласия – статистический критерий для проверки гипотезы о согласии (равенстве) распределения случайной величины исследуемой совокупности с теоретическим распределением или гипотезы о согласии распределений в двух и более совокупностях.

Кстати, любой статистический критерий является случайной величиной, поэтому статистические гипотезы носят вероятностный характер. Уровень значимости \alpha – это вероятность отвержения верной нулевой гипотезы.

\beta – вероятность принятия неверной нулевой гипотезы. Мощность критерия 1-\beta – вероятность отвержения нулевой гипотезы, когда она неверна.

Часто используются следующие критерии согласия:

  1. Критерий Пирсона, или критерий хи-квадрат (\chi^2) (англ. chi-square test, который по сути вычисляет меру различия между гистограммой полученных данных и гистограммой теоретического распределения, а далее сравнивается полученное значение с предполагаемым значением (квантили можно взять из таблицы или посчитать самим) при заданных уровне значимости и числе степеней свободы df. Если полученное число меньше предполагаемого, то распределение подходит.

    df вычисляется как число интервалов минус число оцениваемых параметров распределения и минус 1. Например, если сравниваем с нормальным распределением, то число оцениваемых параметров равно двум (дисперсия и среднее), значит df=n-2-1=n-3

    \chi^2* = n*\sum_{i=1}^{m}\frac{(p_i*-p_i)^2}{p_i}

    Есть соответствующие функции в scipy, matlab.

  2. Критерий Колмогорова-Смирнова (англ. KS-test) – похож на критерий Пирсона, только анализируется разность между функцией распределения (cdf) предполагаемого распределения и имеющегося распределения данных измерений.

    Есть соответствующие функции в scipy, matlab.

Попробуем применить критерий согласия к нашим данным.

Примем за нулевую гипотезу распределение Эрланга. Проверим эту гипотезу используя критерий согласия Пирсона. В коде ниже происходит следующее:

  1. В переменную p_observed записываем частоты появления интервалов. В переменную bins – границы интервалов.
  2. В tmp записываем значения функции распределения Эрланга для наших интервалов. В качестве параметров распределения берем значения из предыдущего шага с графическим методом. Тут я немного сместил центр распределения (параметр loc) для более точного подбора распределения.
  3. В цикле считаем вероятности значений в рамках наших интервалов. Далее получаем частоты появления интервалов согласно предполагаемому распределению Эрланга
  4. С помощью функции scipy получаем вероятность того, что нулевая гипотеза верна. Получилось около 10%.
  5. Далее делаем всё тоже самое, что делает функция scipy, только сами.

В прилагаемом файле ipython notebook вы также найдете некоторые комментарии.

counts, bins = np.histogram(latencies, bins=bins_number)
p_observed = np.array(counts, dtype=float)
#print(p_observed)

tmp = ss.erlang.cdf(bins, k, loc=erlang_loc-0.55, scale=lambd)
#print(tmp)

# Now get probabilities for intervals
p_expected = []
for idx, p in enumerate(tmp):
    if idx == len(tmp) - 1: break
    p_expected.append(tmp[idx+1] - tmp[idx])

p_expected = n*np.array(p_expected)

# Calculate using scipy built in chi-square test
chi_test_value, probability_fit = ss.chisquare(p_observed, np.array(p_expected, dtype=float))
print('Вероятность того, что распределение Эрланга подходит к полученным данным: {0:f}'.format(probability_fit))

# Calculate by ourselved
chi_star = np.sum((p_observed-p_expected)**2/p_expected)
print("chi_star = {}".format(chi_star))
conf_interval = 0.95
df = len(bins) - 3
chi = ss.chi2.ppf(conf_interval, df)  # chi-square quntile for alpha = conf-interval, degrees of freedom = df
print("chi = {}".format(chi))

    Вероятность того, что распределение Эрланга подходит к полученным данным: 0.098665
    chi_star = 27.2634383493
    chi = 28.8692994304

\chi>\chi* – это указывает на то, что распределение Эрланга подходит для наших данных.

Однако, стоит отметить, что разница между величинами небольшая.

В коде также есть вычисление квантиля распределения (функция ppf, percent point function). Квантиль – это значение обратной функции распределения. Если в обычной функции распределения F(x) в качестве x мы принимаем какое-либо значение и получаем вероятность того, что случайная величина X с этим распределением примет значения от минус бесконечности до x, то в обратной функции мы решаем обратную задачу: при заданном p мы получаем x, то есть получаем такое значение x, которое случайная величина не превысит с заданной вероятностью p.

Если вы предполагаете, что у вас нормальное распределение, вы можете воспользоваться еще одним простым способом проверки – методом моментов. В нём просто сравниваются выборочное среднее значение (среднее значение по полученным данным) с математическим ожиданием предполагаемого распределения и тоже самое делается для дисперсии. Чем меньше разница – тем лучше.

Итог

Сегодня мы так и не добрались до доверительных интервалов, однако изучили много нового: основные определения статистики c распространенными синонимами, два основных подхода выбора подходящего распределения для экспериментальных данных: графический и статистический, а также попробовали всё это на практике, используя превосходных open-source пакет NumPy/SciPy для Python.

Скачать файл ipython notebook

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

Используемая литература

  1. http://www.statsoft.com/textbook/distribution-fitting/ – много занимательных gif-ок, в которых показывают различные распределения и измерения кривых распределения и плотности в зависимости от различных значений параметров распределений. Полезная информация для графического подбора распределения.
  2. http://iglin.exponenta.ru/All/ContData/ContData.html – полезная информация о проверки стат. гипотез.
  3. http://sernam.ru/book_tp.php?id=19
  4. http://vfkomd.ru/docs/lections/ms/18_ocenka%20raspred.pdf – глава из неизвестного учебника про подбор распределения.
  5. http://www.hr-portal.ru/statistica/gl3/gl3.php – просто и понятно про статистику
  6. http://www.mathworks.com/... – немало информации по теме от mathworks.
  7. http://www.slideshare.net/enthought/numpyscipy-statistics – отличная презентация о возможностях nympy/scipy.
  8. http://www.grandars.ru/student/statistika/generalnaya-sovokupnost.html – хорошее практическое объяснение всех терминов.
  9. http://statanaliz.info/teoriya-i-praktika/10-variatsiya/15-dispersiya-standartnoe-otklonenie-koeffitsient-variatsii.html – статистика простыми словами
  10. http://esotericpl.narod.ru/Ps_enc/Ps_enc9_4_0.html – хорошо рассказано про статистику малых выборок
  11. http://www.enviroliteracy.org/pdf/materials/1210.pdf – хорошая работа по критерию Пирсона.