Архив рубрики: Python

Исходный код eco.asmon.ru

Выкладывают в свободный доступ исходный код проекта Карты сдачи вторсырья о котором я как-то писал. Код может быть интересен как пример использования Google Fusion Tables, Google Maps и Django.
Код: https://github.com/DimitryDushkin/eco_map.

Интерактивная карта пунктов приёма вторсырья.

Карта пунктов приема вторсырья

В нашей лаборатории сделали небольшой проект – интерактивную карту пунктов приёма вторсырья – Переработай это! –  eco.asmon.ru.

Интерактивная карта пунктов приёма вторсырья.

Под капотом bootstrap, django, Google Maps, в качестве БД используется Google Fusion Tables. БД доступна всем для просмотра. Через приложение можно добавить новые точки приёма вторсырья.

В планах:

  1. Адаптация веб-приложения для планшетов и смартфонов;
  2. Геймификация!)

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

Matplotlib: самое короткое решение проблемы отображения кириллицы

В супер-библиотеке рисования графиков Matplotlib есть досадный баг – кириллические подписи в графике отображаются некорректно. В интернете немало решений, но все они какие-то длинные и не pythonic-way.) Вот самое короткое:

from matplotlib import rc, plot

font = {'family': 'Verdana',
        'weight': 'normal'}
rc('font', **font)

raw_data_filename = 'data_classify/raw_data_1.csv'
X_raw = np.loadtxt(raw_data_filename, delimiter=",")

xlabel(u'Скорость 1')
ylabel(u'Скорость 2')
title(u'Значения показателей')
plot(X_raw[:, 0], X_raw[:, 2], '+')
grid()

savefig('test_1.eps')

Получим что-то типа такого:

Все и в eps корректно сохранится, и покажет правильно.

Erlang: аналогии для понимания циклов

Пришел в голову интересный способ объяснения циклов на Erlang:

Пусть А – массив:

A = [1, 2, 3, 4, 5].

Надо пройтись по каждому элементу массива и увеличить его на 1. На обычных языках программирования это делается примерно так (python):

for i in range(len(A)):
   A[i] += 1

PHP

foreach ($A as $i=>$value)
    $A[i]++;

А теперь эквивалент на Erlang:

A2 = lists:map(
   fun(Elem) ->
      Elem + 1
   end, A).

На Erlang мы использовали анонимную функцию, которая в качестве аргумента принимает одно значение и возвращает его, увеличенным на 1. Также мы ввели новый массив A2, т.к. на Erlang нельзя переопределять единожды присвоенное значение переменной. Да, по сути в Erlang переменных нет =), потому и нет классических циклов. Поэтому для реализации циклов используют или рекурсию, или проход по элементам списка с помощью функций lists:map/2, lists:foldl/3 и др.

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

Доверительный интервал

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

Оценка какого-либо параметра распределения с помощью доверительного интервала называют интервальной оценкой.

Главное преимущество интервальной оценки – возможность характеризовать уверенность в вычисленных параметрах распределения. Например, можно будет сказать: «С вероятностью 95% среднее значение \mu нормального распределения для моих данных будет лежать в интервале от 100 до 102, или, как еще можно это записать, \mu=100\pm1.

Зачем вообще нужна интервальная оценка? Предположим, что мы хотим узнать время обработки запроса каким-либо веб-сервисом. Для этого мы производим несколько измерений: мы последовательно посылаем запросы, получаем ответы на запросы, записываем время между отправкой запроса и получением ответа (время отлика). Как вы уже могли заметить из предыдущих записей, в ходе измерения мы получаем множество различных времен отклика. Какие-то значения встречаются чаще, какие-то – реже. Всё это говорит о том, что мы имеем дело со случайной величиной. У случайной величины есть некий закон распределения, который задается рядом параметров. Параметры распределения мы не знаем, но можем их оценить на основе имеющейся выборки (результатов измерений). Например, оценкой дисперсии нормального распределения будет выборочная дисперсия.

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

Согласно центральной предельной теореме при большом количестве измерений наша оценка параметра распределения имеет распределение, близкое к нормальному. Поэтому, чем больше у нас будет выборка, тем ближе выборочная дисперсия будет к дисперсии (параметру) нормального закона распределения! Но когда остановится? Понятно, что если мы сделаем 10 000 измерений, то мы достаточно (с уверенностью в 99,9999%) точно сможем оценить истинные значения параметров распределения, но часто такой возможности нет. На помощь приходит интервальная оценка.

Мы можем сказать так: «Мы не хотим делать 10 000 измерений, но нам достаточно быть уверенным в наших оценках параметров распределения на 95%». Проведя вычисления доверительного интервала, вполне может оказаться, что для 95% уверенности достаточно всего лишь 50 измерений. Существенная экономия сил и ресурсов!)

Обозначим за (1-\alpha)*100\% доверительную вероятность, или надежность, или доверительный уровень (англ. level of confidence) – вероятность того, что истинное значение измеряемой окажется в доверительном интервале. Саму величину \alpha называют уровнем значимости. Математически 1-\alpha можно записать так:
P(a\le \theta \le b)=1-\alpha,
где \theta – это истинное значение оцениваемого интервала, a – нижнее значение интервала, b – верхнее.

Интересно, что \alpha уже встречался ранее, также назывался уровнем значимости, но обозначал вероятность отвержения верной нулевой гипотезы.

Итак, чем больше \alpha, тем меньше доверительная вероятность и тем уже доверительный интервал. Обычно достаточно \alpha=0.05\%, тогда доверительная вероятность будет равна 95%.

У разных распределений разные параметры и, соответственно, построение интервальных оценок для этих параметров разное.

Доверительные интервалы нормального распределения

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

\bar{X}\pm Z_{1-\alpha}\frac{\sigma}{\sqrt{n}}, где n – объем выборки, \overline{X} – выборочное среднее и, самое интересное, Z_{1-\alpha} – квантиль стандартного нормального распределения (нормальное распределение с мат. ожиданием = 0 и стандартном отклонением = 1) для заданной вероятности 1-\alpha.

Не будем приводить все выкладки, почему именно так определяется доверительный интервал. Подробнее вы можете почитать в [1] (стр. 50-55).

Поясню еще раз про квантиль. Квантиль – это такое значение x, которое случайная величина не превысит с заданной вероятностью 1-\alpha. Квантиль – это значение обратной функции распределения, т.е. 1-\alpha=F(Z).

Статистика малых выборок

В случаях, когда стандартное отклонение неизвестно, а также объем выборки не превышает 10-30 значений, используется t-статистика для оценки значения параметра распределения. t-статистика вычисляется как квантиль t-распределения (распределение Стьюдента, см. первую часть этой серии записей) и задается только одним параметром df – степенью свободы. Для нормального распределения df=n-1, где n – объем выборки.

Тогда интервальная оценка мат. ожидания будет выглядеть следующим образом:

\bar{X}\pm t_{\alpha,df} \frac{S_x^2}{\sqrt{n}}, где t_{\alpha,df} – коэффициента Стьюдента (значение квантили статистики t порядка P=1-\alpha/2 для числа степеней свободы df = n – 1), S_x^2 – выборочная дисперсия.

Доверительные интервалы экспоненциального распределения

После предыдущей записи, я попробовал еще одно распределение для имеющихся данных – экспоненциальное. Оказалось, что оно лучше подходит, чем распределение Эрланга. Результаты измерений вы найдете на иллюстрации ниже и в прилагаемом файле ipython notebook. Видимо, это произошло из-за того, что распределение этих конкретных данных не имеет левого хвоста.

Поэтому, продолжая наш пример, найдем доверительный интервал для параметра \lambda экспоненциального распределения. Кстати, если бы мы использовали экспоненциальное распределение в расчетах в рамках теории массового обслуживания, то параметр \lambda называли «интенсивностью поступления заявок» и он обозначал бы среднее время поступления новых заявок за единицу времени.

Оценкой \lambda является обратное значение среднего арифметического по выборке, т.е.
\hat\lambda=1/\bar{x}.

Доверительный интервал выглядит следующим образом:
$latex \frac{2n}{\widehat{\lambda}\chi^2_{1-\alpha/2,2n}} < \frac{1}{\lambda} < \frac{2n}{\widehat{\lambda}\chi^2_{\alpha/2,2n}}$, где $latex \chi^2_{\alpha/2,2n}$ – квантиль распределения хи-квадрат при заданных доверительной вероятностью и степенью свободы $latex k=2n$.

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

#=========================================
# 1) Графический метод
#=========================================

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

print(ss.expon.fit(latencies))
exp_loc, exp_scale = ss.expon.fit(latencies)        # Оценим положение распределения

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

# Получим соответствующие значение по функции плотности вероятности для распределения Эрланга
# В качестве начала построения (параметр loc) выберем минимальное значение  из выборки,
# т.е. сместим распределение так, чтобы оно начиналось с нуля.
# В scipy параметр scale = 1 / lambda, где lambda - параметр экспоненциального распределения.
distr_data = ss.expon.pdf(x, loc=np.min(latencies), scale=4.8)

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


#=========================================
# 2) Точечная оценка параметра распределения
#=========================================


# Оценим с помощью хи-критерия соответствие распределения
counts, bins = np.histogram(latencies, bins=bins_number)
p_observed = np.array(counts, dtype=float)

tmp = ss.expon.cdf(bins, loc=exp_loc-0.6, scale=4.8)
#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 = number_of_samples*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))

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

Следующий код проводит интервальную оценку параметра:

#==================================================
# 3) Интервальная оценка параметра распределения
#==================================================


n = size(latencies)
alpha = 0.05

# Сдвинем данные так, чтобы они начинались с 0
data = np.array(latencies) - np.min(latencies)

# Посчитаем оценку лямбда
lambda_estimate = 1/np.mean(data)
print("Оценка 1/lambda на основе выборки: {}".format(1/lambda_estimate))

# Посчитаем левую границу интервала
a = 2 * n / (lambda_estimate * ss.chi2.ppf(1 - alpha / 2, 2 * n))

# Правую
b = 2 * n / (lambda_estimate * ss.chi2.ppf(alpha / 2, 2 * n))
             
print("95% доверительный интервал: {} < 1/lambda < {}".format(a, b))


ci_width = b - a
print("Ширина доверительного интервала: {}".format(ci_width))

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

    Оценка lambda на основе выборки: 0.249376558603
    Оценка 1/lambda на основе выборки: 4.01
    95% доверительный интервал: 3.32700158323 < 1/lambda < 4.92847012339

Я везде привожу \frac{1}{\lambda}, т.к. это значение использовалось в графическом методе вместо просто \lambda.

Определения необходимого объема выборки для нормального распределения при известном стандартном отклонении

А теперь разберемся с самым интересным! Вот ради чего всё и затевалось.)

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

Опять начнем с простой ситуации, когда известно стандартное отклонение. Такой вариант маловероятен, но тем не менее возможен. Вспомним запись мат. ожидания с доверительным интервалом:
\bar{X}\pm Z_{1-\alpha} \frac{\sigma}{\sqrt{n}}
Величина, добавляемая и вычитаемая из \bar{X}, называется ошибкой выборочного исследования и обозначается следующим образом:
e=Z_{1-\alpha} \frac{\sigma}{\sqrt{n}}

Из формулы e можно определить объем выборки n:
n=\frac{Z^2 \sigma^2}{e^2}

Таким образом, для определения объема выборки необходимо знать три параметра.
1. Требуемый доверительный уровень, который влияет на величину Z, являющуюся критическим значением стандартизованного нормального распределения.
2. Приемлемую ошибку выборочного исследования e. Задается в измеряемых единицах. Например, если бы у нас было нормальное распределение e можно было бы взять равным 5 (мс).
3. Стандартное отклонение ?.

Определения необходимого объема выборки для нормального распределения при неизвестном стандартном отклонении

Здесь, как и при определении доверительного интервала, воспользуемся t-распределением и выборочной дисперсией S_x^2.

n\ge t^2_{\alpha,df}  (\frac{S_x}{\delta})^2=t^2_{\alpha,df} \epsilon^2,
где \delta – требуемая точность. Обычно требуемую точность берут в размерах выборочной дисперсии. Например, если нам нужна точность измерений в 2 раза выше выборочной дисперсии, то \epsilon=\frac{S_x}{0.5S_x}=2.

Стоит отметить, что по сути выборочная дисперсия зависит от объема выборки, поэтому вычисляемый объем выборки необходимо решать методом последовательных приближений. В качестве первого значения выборки можно взять, например, df=10. Посчитать n. Получим, например, 8. Снова n, но при df=8-1=7. Получим n=7. И т.д., пока значение n станет приблизительно неизменным.

Определения необходимого объема выборки для экспоненциального распределения

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

Поставим следующую задачу: сколько необходимо измерений, чтобы быть уверенным в том, что оценка параметра \lambda попадет в доверительный интервал с вероятностью 95%?

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

Из неравенства доверительного интервала выведем формулу расчета необходимого объема выборки. Значение ширины доверительного интервала рассчитываем как разницу между верхней и нижней интервальными оценками и примем за \epsilon (эпсилон). Это значение мы зафиксируем.

\frac{1}{n} = \frac{2}{\epsilon \widehat{\lambda}} (\frac{1}{\chi^2_{\alpha/2,2n}} - \frac{1}{\chi^2_{1-\alpha/2,2n}})

Код ниже просто рассчитывает приведенную формулу:

#==================================================
# 4) Определение необходимого объема выборки 
# при фиксированной ширине доверительного интервала
#==================================================

n = 100
epsilon = ci_width
alpha = 0.1

n_ci = (2 / (epsilon * lambda_estimate)) * (1/ss.chi2.ppf(alpha / 2, 2 * n) - 1/ss.chi2.ppf(1 - alpha / 2, 2 * n))
n_ci = 1 / n_ci
print('Необходимый объем выборки = {} для попадания истинного значения параметра распределения' \
     ' в интервал шириной = {} с доверительной вероятностью = {}%'.format(round(n_ci), epsilon, (1-alpha)*100))

Результ: «Необходимый объем выборки = 120.0 для попадания истинного значения параметра в интервал шириной = 1.60146854016 с доверительной вероятностью = 90.0%».

Нахождение объема выборки при фиксированной доверительной вероятности и ширине доверительного интервала

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

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

Код (который вы также найдете в файле ipython notebook в конце поста):

#==================================================
# 5*) Определение необходимого объема выборки 
# при фиксированной ширине доверительного интервала и доверительной вероятности
#==================================================
alpha = 0.05
epsilon =  0.5 * ci_width
optimum_n = 0
N = [0]

for i in range(0, 500, 10):
    n = i
    
    for k in range(5):
        n_ci = (2 / (epsilon * lambda_estimate)) * (1/ss.chi2.ppf(alpha / 2, 2 * n) - 1/ss.chi2.ppf(1 - alpha / 2, 2 * n))
        n_ci = 1 / n_ci
        n = n_ci
        N.append(n_ci)
        #print("{}".format(n_ci))
    
    #print('---')
    
    # Замечено, что при увеличении первоначального значения n
    # в определенной итерации метода приближения вычисленный объем выборки
    # перестает падать или расти, а колеблется около одного и того же числа.
    # Это число можно считать оптимальным объемом выборки при заданной ширине доверительного интервала.
    # Что требует доказательства.
    
    # Ниже мы грубо оцениваем размах колебания значения n_ci в этой итерации.
    # Если оно мало, значит это значение n можно считать оптимальным.
    if np.std(N[-5:]) < 5:
        optimum_n = n_ci

print('Эверстически оптимальный объем выборки: {}'.format(math.ceil(optimum_n)))
print('При заданной ширине доверительного интервала: {}'.format(epsilon))
print('И доверительной вероятности: {}%'.format((1-alpha)*100))
plot(N)

Если использовать в качестве начального значения n изначальный объем выборки (в нашем случае это 100), то каждое следующее приближение будет уменьшать вычисленной значение n_ci до нуля, однако, я заметил, что при увеличении первоначального значения n (больше первоначального объем) в определенной итерации метода приближения вычисленный объем выборки перестает падать или расти, а колеблется около одного и того же числа. Это число можно считать оптимальным объемом выборки при заданной ширине доверительного интервала. Это, конечно, требует доказательства и является полнейшей эвристикой.

Полученный с помощью этой эвристики результаты хорошо согласуются с предполагаемым поведением этой функции: чем шире значение доверительного интервала, тем меньше нужен объем выборки и наоборот. Однако сам алгоритм нестабилен и не выдает результатов при \epsilon на порядок меньше исходной ширины интервала. При больших значениях он соответственнно выдаст 0 (что можно понимать как выборку с 1 значением), т.к. если интервал достаточно широк, то попасть в него будет очень просто с любой выборкой.

Ниже график, на котором показано вычисление объема выборки методом последовательных приближений при различных начальных значениях n. Хорошо виден интервал, в котором вычисленное значение колеблется, но заметно не растет и не падает.

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

Скачать файл с ipython notebook раздела

Заключение

Занимаясь написанием этой череды постов, я узнал много нового про распределения и методы работы с ними, надеюсь и вы почерпнули для себя что-нибудь интересное и полезное. =)

Список литературы

  1. Спирин H.A. Методы планирования и обработки результатов инженерного эксперимента : Конспект лекций. / H.A. Спирин, В.В. Лавров. Екатеринбург : ГОУ ВПО УГТУ-УПИ, 2004. — 257 с.
  2. Статистика для менеджеров с использованием Microsoft Office Excel 4-e издание. Дэвид М. Левин и др. – отлично рассказано в 7 главе про практическое использование доверительных интервалов. И именно 7-я глава есть в открытом доступе.
  3. Confidence Intervals, Guy Lebanon – отличная работа, рассказывающая о доверительных интервалах.
  4. http://en.wikipedia.org/wiki/Exponential_distribution#Maximum_likelihood – про оценку параметра экспоненциального распределения.

Разберемся с измерением случайных величин (ч. 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 – хорошая работа по критерию Пирсона.

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

Основная терминология

Поскольку теория вероятности нашла применение во многих областях, появилось немало обозначений для одного и того же. Ниже я выпишу русские и английский термины с их синонимами. Синонимы в рамках одного языка пишутся через запятую. Знание английского варианта может быть очень полезно при работе с различными математическими системами (типа Matlab, Mathematica) и библиотеками.

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

Distribution function, Сumulative distribution function, c.d.f.  Функция распределения, Функция распределения вероятностей, Интегральная функция распределения вероятностей – функция, характеризующая распределение случайной величины. Записывается как F(x) и обозначает вероятность того, что после выполнения эксперимента случайная величина X получит значение, не превышающее число x. Математически это записывается как
F(x) = P(X \leq x), -\infty < x < +\infty

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

Дискретная (от лат. disctretus – прерывистый) случная величина может принимать только определенный набор значений.

Например, результат бросания монеты может быть только орел или решка.

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

Возвращаясь к примеру из первой части, время открытия веб-страницы – это непрерывная случайная величина.

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

Probability massProbability mass functionp.m.f. – Функция вероятности – функция, описывающая распределение дискретной случайной величины. Записывается как множество значений, которые может принять случайная величина, и соответствующие им вероятности. Например, для монеты p_1=0.5, p_2=0.5, где p_1 – вероятности выпадения орла и p_2 – решки.

С функцией вероятности у меня остается один непонятный момент. Почему функцию распределения плотности вероятности для дискретной случайной величины назвали именно «функция вероятности»? Это же очень похоже на изначальное понятие «функция распределения вероятностей»? Надеюсь, кто-нибудь из читателей меня просветит.

Впрочем, дальше будет речь идти только о непрерывных величинах, поэтому этот вопрос некритичный.

Probability densityProbability density functionp.d.f  – Плотность вероятности, Функция плотности вероятности – функция, описывающая распределение непрерывной случайной величины. Кстати, да, «density» дословно переводится как «плотность». Плотность вероятности записывается как f(x) . Математически плотность вероятности является производной от функции распределения вероятностей:

f(x)=F'(x), -\infty < x < +\infty.

Соответственно, функцию распределения вероятностей можно выразить через плотность:

F(x)=P(X\in[\infty,x])=\int\limits_{-\infty}^{x} f(x)dx

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

Само обозначение «плотность» взято из механики. В этой интерпретации функция f(x) буквально характеризует плотность распределения масс по оси абсцисс (ось x):

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

Виды распределений

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

  1. F(x)\in[0;1], потому что не может быть отрицательной вероятности и вероятности больше единицы;
  2. F(x_1)\le F(x_2),x_1<x_2, т.е. функция неубывающая, потому что чем больше совокупность случайных событий, тем более вероятно их наступление;
  3. \lim_{x\to -\infty}=0, \lim_{x\to +\infty}=1 по той же причине, что и п.2.
В википедии можно найти по меньшей мере 30 различных законов распределения случайных величин. Разберем одно из самых популярных распределений – нормальное, или распределение Гаусса.

Нормальное распределение

О нем говорилось в прошлой части, но теперь, когда вам знакомы понятия функции и плотности распределения, мы можем рассмотреть его в новом свете!
Плотность распределения записывается следующим образом:
В принципе, саму формулу запоминать нет особого смысла, главное запомнить, что в ней присутствуют два параметра, которые влияют на форму распределения (см. график ниже): \mu (мю) – среднее значение (математическое ожидание), указывает «центр» «купола», и \sigma^2 (сигма) – дисперсия (рассеивание), влияет на «толщину» купола. Сигму не в квадрате называют среднеквадратическим отклонением.
Выглядит плотность вероятности f(x) так:
Рассматривая на графике красную линию, можно сказать, что наиболее вероятное значение – это ноль, вероятность которого примерно равна 0.9.
Используя интерактивную консоль ipython с поддержкой построения графиков мы можем легко построить такие же графики сами (надеюсь, вы установили дистрибутив The Enthought Python Distribution из предыдущей части)!
Для начала запустите ipython в консоли командой «ipython —pylab». Затем введите следующий код:
from scipy.stats import norm
import numpy as np

x = np.linspace(-5,5,100)
plot(x, norm.pdf(x))
Вы получите следующий график:

Краткий итог

Да, эта часть получилась теоретическая… Зато она вносит ясность в самые базовые элементы статистики и теории вероятности, которые обязательно понадобятся нам в будущем! А для практикующих, надеюсь, была полезна часть с разбором синонимов в терминологии.

Смотрите в следующей серии

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

ipython + pylab (matplotlib) в составe Enthought Python Distribution

Проблема: русские буквы отображаются прямоугольниками на графиках matplotlib.

Решение: необходимо задать подходящий backend для matplotlib. Для этого:

  1. Узнаем, где находится файл конфигурации matplotlib:
    >>> import matplotlib
    >>> matplotlib.matplotlib_fname()
    '/home/foo/.matplotlib/matplotlibrc'
    
  2. Находим в нем строку «backend: …»
  3. Если у вас MacOS заменяем на «backend: MacOSX», если что-то другое, то попробуйте что-нибудь другое, например, GTKAgg или Qt4Agg. Полный список можно найти здесь.