Архив метки: python

Наука: наш доклад на АКТО-2012

Вчера были в Казани – пожалуй, единственном городе в России (из тех, что я видел), где архитекторам не ограничивают фантазию!)

Впрочем, мы были там по делу – выступали на конференции «Авиакосмические технологии, современные материалы и оборудование. Казань-2012» в секции 7 – «Современные требования к подготовке кадров для работы в авиакосмической отрасли. Новые образовательные технологии в подготовке и переподготовке кадров.»

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

Основная заслуга по праву принадлежит Наталии Смирновой – главному теоретическому и практическому исследователю и разработчику системы.) Я в основном консультирую и помогаю по технической стороне вопроса, именно поэтому проект разрабатывается на Django 😉

Разберемся с измерением случайных величин (ч. 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%!)

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

Задача: определить, сколько раз необходимо провести измерения производительности веб-сервиса, чтобы быть уверенным в полученных результатах на 95%.

Эта серия записей не претендует на звание статьи, методического пособия или еще чего-то. Я хочу рассказать об измерениях случайных величин просто, «на пальцах». Так, чтобы даже без знания теории вероятности после прочтения записи было бы всё понятно и можно было бы попробовать свои силы на практике.

Для начала необходимо признать один простой факт – результат любого измерения – это случайная величина. В моём случае результат измерения среднего времени обработки запросов веб-сервисом – случайнее вдвойне!)

Тут в дело не только в погрешности оборудования клиента и сервера, сети, но и во времени суток! В «час пик» среднее время обработки запроса может быть значительно больше. И тут в игру вступают условные вероятности. Впрочем, стоп. Сегодня пройдемся только по простому измерению без всяких условных факторов.

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

# coding:utf-8

import httplib
import time

for i in range(10):
    conn = httplib.HTTPConnection("ag.ru")
    conn.request("GET", "/")
    start_time = time.time()
    r1 = conn.getresponse()
    latency = time.time() - start_time
    print("Время обработки запроса: {0}мс".format(round(latency * 1000)))

И получите примерно такой результат:

Время обработки запроса: 14.0мс
Время обработки запроса: 10.0мс
Время обработки запроса: 10.0мс
Время обработки запроса: 10.0мс
Время обработки запроса: 15.0мс
Время обработки запроса: 9.0мс
Время обработки запроса: 14.0мс
Время обработки запроса: 9.0мс
Время обработки запроса: 10.0мс
Время обработки запроса: 11.0мс

Обозначим прописной буквой X возможное время обработки запроса, маленькой же буквой x обозначим конкретные принимаемые значения времени обработки запроса. Для приведенного примера для сайта http://ag.ru X может принимать значения от 9 до 15 мс, для первого измерения x_1=14мс.

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

Для запуска нового кода вам понадобится установленная библиотека matplotlib, которую в простейшем случае можно установить с помощью команды в консоли pip install matplotlib. Я же рекомендую скачать полный набор необходимых библиотек, который содержится в Enthought Python Distribution (тут и numpy, и matplotlib, и еще немало полезных и часто используемых библиотек), дистрибутив которого разработчики подготовили для всех популярных платформы (Windows, MacOS, Linux).

Код обновленной программы будет следующий:

# coding:utf-8

import httplib
import time
import matplotlib.pyplot as plt

latencies = []
for i in range(100):
    conn = httplib.HTTPConnection("habrahabr.ru")
    conn.request("GET", "/")
    start_time = time.time()
    r1 = conn.getresponse()
    latency = round((time.time() - start_time) * 1000)
    latencies.append(latency)
    print("Время обработки запроса: {0}мс".format(latency))

plt.hist(latencies, 50)
plt.title(u"Гистограмма времен обработки запросов")
plt.xlabel(u"Задержка")
plt.ylabel(u"Частота")
plt.show()

И получим следующий график:

Гистограмма времен обработки запросов для сайта habrahabr.ru

Какие же выводы мы можем сделать из этого графика? Получается, чаще всего вы сможете увидеть главную страницу сайта habrahabr.ru где-то спустя 900 мс, вероятность такого события самая высокая, поскольку именно за 900 мс +- несколько десятков миллисекунд согласно нашим измерениям мы 20 раз получали главную страницу! Также есть небольшие вероятности, что страница загрузится или за 100 мс, 200 мс и т.д.

Наша гистограмма весьма похожа на нормальное распределение! Нормальное распределение – это такое распределение случайной величины, где есть единственное наиболее вероятное значение случайного события – среднее значение, или математическое ожидание. Вероятности принятия случайной величиной других значений плавно убывают симметрично относительно среднего значения. Нормальным его называют потому, что именно такое распределение чаще всего встречается в природе. Например, рост человека – это случайная величина с нормальным распределением. Чаще всего вы встретите на улице человека с ростом 170-180 см, но есть же и люди с ростом 155 см и под два метра.

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

Конечно, разные веб-сервисы будут по разному выглядеть на такой гистограмме, и здесь нам весьма повезло, что получили близкую к канонической форме «купола» нормального распределения.

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

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

  1. Кельтон В., Лоу А. Имитационное моделирование. Классика CS. 3-е изд. – СПб.: Питер; Киев: Издательская группа BHV, 2004.-847 c:ил.
  2. Спирин H.A. Методы планирования и обработки результатов инженерного эксперимента : Конспект лекций. / H.A. Спирин, В.В. Лавров. Екатеринбург : ГОУ ВПО УГТУ-УПИ, 2004. — 257 с.

Используемое программное обеспечение (ПО)

  1. Язык программирования Python
  2. Библиотека построения графиков matplotlib в составе дистрибутива The Enthought Python Distribution

Google API Python Client: сохранение access token при использовании service account

Проблема: текущий python клиент Google API не хранит access token в режиме доступа к API с помощью service account, что происходит из-за отсутствия полноценной поддержки хранилища данных (oauth2client.file.Storage). При каждом новом запросе к функциям API необходим дополнительный запрос на обновление access token, что занимает немало времени (около секунды).

Решение: изменение исходного кода oauth2client.client.

Теперь следующий код запроса к функциям API будет работать:

import httplib2
from apiclient.discovery import build
from oauth2client.client import SignedJwtAssertionCredentials
from oauth2client.file import Storage

logging.getLogger().setLevel(logging.INFO)

f = file('616591315bc79b4b964cf6c2d4afbe3b1a82a1c1-privatekey.p12', 'rb')
key = f.read()
f.close()

# First argument is e-mail provided in Google API Console for service account
credentials = SignedJwtAssertionCredentials(
    '...@developer.gserviceaccount.com',
    key,
    scope='https://www.googleapis.com/auth/fusiontables')
storage = Storage('fusion.dat')
credentials.set_store(storage)

http = httplib2.Http()
http = credentials.authorize(http)

service = build("fusiontables", "v1", http=http)
response = service.query().sqlGet(sql='SELECT Location FROM 1gvB3SedL89vG5r1128nUN5ICyyw7Wio5g1w1mbk LIMIT 10').execute(http)
for row in response['rows']:
    print unicode(row[0]).encode('utf-8')

При текущей версии клиента код выдавал ошибку:

UnicodeDecodeError: 'utf8' codec can't decode byte 0x82 in position 1: invalid start byte

Я немного поковырялся в исходниках и исправил ситуацию. Исходный файл пропатченного client.py для пакета oauth2client прикрепляю к этой записи. Коммит пока не делал, т.к. не до конца уверен, что изменения не затронут другие режимы (client, web-application).

PS: кстати, если у вас при использовании клиента даже в простом режиме появляется ошибка «invalid_grant» – проверьте системное время! Я долго бился об стенку, пока на нашел такое предположение, и оно оказалось верным! Может быть, это из-за second leap?! Не скажу точно.

Скачать пропатченный исходный код oauth2client.client.

Google API Python Client: Google Fusion Tables

Наконец-то я получил доступ к google fusion tables API с помощью google-api-python-client. У этой библиотеки действительно странноватый интерфейс, но всё стало гораздо понятней после недавнего полноценного переезда Fusion API на новую версию.

Рабочий код:

import httplib2
from apiclient.discovery import build
from oauth2client.client import SignedJwtAssertionCredentials

# Here's the file you get from API Console -> Service Account.
f = file('key.p12', 'rb')
key = f.read()
f.close()

# Create an httplib2.Http object to handle our HTTP requests and authorize it
  # with the Credentials. Note that the first parameter, service_account_name,
  # is the Email address created for the Service account. It must be the email
  # address associated with the key that was created.
credentials = SignedJwtAssertionCredentials(
    '...@developer.gserviceaccount.com',
    key,
    scope='https://www.googleapis.com/auth/fusiontables')
http = httplib2.Http()
http = credentials.authorize(http)

service = build("fusiontables", "v1", http=http)
# For example, let make SQL query to SELECT ALL from Table with
# id = 1gvB3SedL89vG5r1128nUN5ICyyw7Wio5g1w1mbk
print(service.query().sql(sql='SELECT * FROM 1gvB3SedL89vG5r1128nUN5ICyyw7Wio5g1w1mbk').execute())

Можно использовать функцию help, чтобы получить документацию по методам:

service = build("fusiontables", "v1")
help(service)
help(service.query())

Для последней команды вы получите что-то типа такого:


class Resource(__builtin__.object)
 |  A class for interacting with a resource.
 |
 |  Methods defined here:
 |
 |  __init__(self)
 |
 |  sql = method(self, **kwargs)
 |      Executes an SQL SELECT/INSERT/UPDATE/DELETE/SHOW/DESCRIBE statement.
 |
 |      Args:
 |        pp: string, A parameter
 |        hdrs: boolean, Should column names be included (in the first row)?. Default is true.
 |        typed: boolean, Should typed values be returned in the (JSON) response -- numbers for numeric values and parsed geometries for KML values? Default is true.
 |        strict: string, A parameter
 |        userip: string, A parameter
 |        sql: string, An SQL SELECT/SHOW/DESCRIBE/INSERT/UPDATE/DELETE statement. (required)
 |        trace: string, A parameter
 |
 |      Returns:
 |        An object of the form:
 |
 |          { # Represents a response to an sql statement.
 |          "kind": "fusiontables#sqlresponse", # Type name: a template for an individual table.
 |          "rows": [ # The rows in the table. For each cell we print out whatever cell value (e.g., numeric, string) exists. Thus it is important that each cell contains only one value.
 |            [
 |              [
 |                "",
 |              ],
 |            ],
 |          ],
 |          "columns": [ # Columns in the table.
 |            "A String",
 |          ],
 |        }
 |
 |  sqlGet = method(self, **kwargs)
 |      Executes an SQL SELECT/SHOW/DESCRIBE statement.
 |
 |      Args:
 |        pp: string, A parameter
 |        hdrs: boolean, Should column names be included (in the first row)?. Default is true.
 |        typed: boolean, Should typed values be returned in the (JSON) response -- numbers for numeric values and parsed geometries for KML values? Default is true.
 |        strict: string, A parameter
 |        userip: string, A parameter
 |        sql: string, An SQL SELECT/SHOW/DESCRIBE statement. (required)
 |        trace: string, A parameter
 |
 |      Returns:
 |        An object of the form:
 |
 |          { # Represents a response to an sql statement.
 |          "kind": "fusiontables#sqlresponse", # Type name: a template for an individual table.
 |          "rows": [ # The rows in the table. For each cell we print out whatever cell value (e.g., numeric, string) exists. Thus it is important that each cell contains only one value.
 |            [
 |              [
 |                "",
 |              ],
 |            ],
 |          ],
 |          "columns": [ # Columns in the table.
 |            "A String",
 |          ],
 |        }
 |
 |  ----------------------------------------------------------------------
 |  Data descriptors defined here:
 |
 |  __dict__
 |      dictionary for instance variables (if defined)
 |
 |  __weakref__
 |      list of weak references to the object (if defined)

Up.: нашел список методов доступа к Fusion Tables. В нём вся та же информация, что даёт функция help(), но в виде веб-страниц. Кстати, там есть инфа и по другим API.

Для каждого языка программирования свои задачи, или выводы «9 лет спустя»

Окидывая взглядом последние 9 лет практики программирования, я пришел к ряду интересных и, думаю, полезных выводов. Свою историю изучения языков программирования (далее «ЯП») проиллюстрировал ниже. Где-то отметил фреймворки на этих ЯП. Отмеченные фреймворки, как правило, давали новую мотивацию и открывали новые горизонты в изучении и использовании этих ЯП.

Языки программирования по жизни

Языки программирования по жизни

Вывод 1: для разных классов задач – разные языки программирования

В действительности, достаточно хорошо изучив какой-либо язык, кажется, что на нём надо писать не только профильные приложения. Но будьте осторожны! На этом пути вы можете впасть в состояние «вбивания гвоздей микроскопом». Так рождаются демоны (постоянно висящие процессы, типа серверов) на PHP, веб-сервера на JavaScript, веб-приложения на C++ и прочее. В качестве эксперимента – это прекрасно, но всё становится гораздо ужасней, когда что-то в силу своей популярности в другой области становится частым решением в других.

Хорошим примером иллюстрирующим первый вывод является node.js и Erlang. На node.js с радостью накинулось большое число веб-разработчиков, ведь все знают JavaScript! Попробовав домашние проекты на этой платформе, многие начали использовать node.js в реальных проектах и в итоге столкнулись с большими проблемами, решение которых затрачивало всё больше сил. Например, отказоустойчивость таких приложений. Стоит вылезть ошибке в одном из потоков выполнения, и весь сервер падает, если только вы не оградились везде, где только можно try-catch блоками. Или горизонтальное масштабирование. Тут вообще беда. И рождаются костыли-костыли-костыли.

Другое дело Erlang, который изначально «by-design» создавался для больших телекоммуникационных систем. Писать на нём серверное приложение — сплошное удовольствие. Хотите постоянно работающий сервер, которым можно управлять по telnet? Да легко! Держите встроенный пакет gen_tcp. Горизонтальное масштабирование существующего приложения? Вам скорей всего придется поменять пару строчек, где описывается передача сообщений другим функциям, — добавить в них отправку сообщений на другие ноды. Тоже и с отказоустойчивостью (привет супервизорам), обработкой бинарных данных и пр. пр.

Вывод 2 – не упускайте возможности изучать новые ЯП и/или их популярные фреймворки

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

Существует две основные когорты ЯП:

  1. Си-подобные: С++, С#, Java, PHP, JavaScript, ActionScript и др.
  2. LISP-подобные: Erlang, Python, Ruby, Perl и др.

Если вы знаете хотя бы один язык из каждого класса, то считайте, что вы знаете и все остальные ЯП.) Исключениями можно считать Assembler и языки-приколы, типа Brainfuck. Кстати, часто различить эти классы ЯП можно по наличию/отсутствию кортежей и обрамлению условных операций и циклов в фигурные скобки.

Я не стал разбивать на группы ООП, функциональные и процедурные парадигмы программирования, т.к. по большому счету почти любой подход можно использовать почти в любом языке (например, функциональный подход в Java). Другое дело, что язык изначально для этого не создавался (см. вывод 1;-) и такие извороты выглядят крайне противоестественно. Исключение Python. Он изначально разрабатывался мульти-парадигмальным.

Мой топ-рейтинг открытий типа: «А-а-а! Так вот с помощью чего это надо было делать!»:

1.  Zend Framework (на PHP)

Ну нафига я писал свою глючную библиотеку для работы с БД? О боги, а зачем я каждый раз придумывал новую архитектуру для каждого нового веб-приложения? MVC, Zend_Db, Zend_Action даруют возможность писать веб-приложения быстрее, а код делать надежным!

2.  Erlang

Боже! Ну почему я потратил кучу времени на написание этого простого демона на Java, если в Эрланге это можно было сделать парой десятков строчек, да и работало бы безотказно, да и масштабировать было бы проще!

3.  Octave

Так, надо бы быстрое преобразование Фурье… О, Matlab! Ну что, пошли на страшное преступление в ближайший торрент-трекер. Хотя ладно, вот есть бесплатная библиотечка для любимого PHP. Эхх…что-то как-то всё грустно и медленно. Octave? Octave… Octave! Мега-удобный синтаксис работы с матрицами, встроенные функции для кучи математических задач, да тут еще и графики с пол пинка строятся! Да оно еще совместимо с MatLab! А еще open-source!

4.  KnockoutJS (на JavaScript).

Как же мне надоели эти простыни с кодом, запутанная логика работы с объектами на странице, эти бесконечные onSomeEvent функции. Как бы хотелось MVC, но вот только на клиенте. Чтобы вот есть у тебя модель данных и при её изменении автоматически всё само менялось на странице. Добавил к ней объект «Задача», а на странице возьми и появись соответствующий блок для новой задачи.

KnockoutJS, как ты вовремя!

5.  Django (на Python)

Да-да, Zend Framework, ты не плох, но ты… рутинный. Каждый раз надо писать свою админку, прикручивать авторизацию, свои однотипные операции добавления, редактирования и удаления сущностей, расписывать формы ввода данных, писать однотипные виды для типичных задач, типа показа по 10 новостей на одной странице и разбивках остальных новостей по страницам. Copy-paste, copy-paste…

Ммм… что-то зачастили на хабре со своими дифирамбами во славу Django и Python.  Да давно уже поют, который год. Что ж они всё никак не угомонятся?) Ну попробуем и мы.

Так-так… значит я раскомментирую вот тут и вот тут, еще одна команда в консоли и… у меня готова полноценная админка?! Вы издеваетесь?) Так просто! А вот еще бы хотелось, чтобы в админке было удобно приписывать к разным авторам разные статьи, чтобы было красивая связь многим-ко-многим. Добавить одну строку в модель?! Ну-у, у вас совсем нет совести! Зачем же я столько времени тратил на эту рутину в зенде?»

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