Во всем Хабре сыскалась лишь пара статей на вышеуказанную тему. А тема благодатная. Да и в минувшую среду как раз окончился курс "Introduction to Computational Finance and Financial Econometrics". По мотивам его пятой недели «Descriptive statistics» и появился этот пост. Причастившимся будет неинтересно, а желающих познакомиться с базовыми приемами анализа данных при помощи R — прошу под хабракат.
Предварительные соглашения
О терминах
У автора из статистики был только семестр «тервера» N лет назад. Поэтому после сомнительно переведенных слов и их сочетаний будет указан исходный английский термин (курсивом в скобках). Специалисты, пожалуйста, шлите в личку более корректные варианты терминов. Спасибо.
Об установке
На установке ПО внимание не заостряется намеренно, в виду тривиальности. По крайней мере на Windows платформе все свелось к стандартному «далее -> далее ->… -> готово». Единственный требуемый для выполняемого в статье кода пакет PerformanceAnalytics устанавливается через меню «Пакеты / Установить пакет(ы)...», выбор ближайшего к вам зеркала, выбор нужного пакета из списка.
Набор данных
Хотелось избежать типичности: продажи, квартиры, рентабельность акций (simple returns), — сколько можно? Поэтому предметная область нашей выборки — вечна как в контексте Хабра, так и вне его контекста. Не так давно в блоге СамиЗнаетеКого был опубликован опрос «Какого размера у вас грудь?». Учитывая, что в него были включены два варианта ответа для отсева нерелевантной аудитории, есть некоторая уверенность в правдоподобии выборки. Для удобства результаты приведены и здесь:
Цель
В рамках нашего мини-исследования сравним с нормальным распределением 2 набора данных:
- (НД1) варианты с третьего (уже не выросла) по девятый (шестой размер),
- (НД2) те же варианты, но к нулевому размеру добавим оптимисток (еще не выросла) и включим в данные седьмой размер.
Опытным статистикам очевидно, что изменениями второго варианта мы отдалим распределение от нормального. К финалу статьи у нас накопится достаточно сведений, чтобы обосновать это формально.
Ход исследования
Для начала поместим наши наборы данных в переменные:
data = c(rep(0, 184), rep(1, 510), rep(2, 996), rep(3, 763), rep(4, 327), rep(5, 147), rep(6, 60))
data_ol = c(data, rep(0, 51), rep(7, 65))
x.txt = "Размер груди" # и заодно сохраним повторяющееся наименование оси
Функция c «склеивает» свои аргументы в единый вектор, функция rep(x, y) возвращает вектор из y значений x. Например, rep(0, 184) вернет вектор из ста восьмидесяти четырех нулей. В рекомендациях гугла и еще в нескольких источниках встречалось мнение, что символ равенства негоже использовать для присваивания, лучше — "<-". Знающие люди, пожалуйста, изложите в комментариях достаточно веские обоснования, чтобы писать 2 символа вместо одного. Лично для автора эта альтернатива отдает неудобством оператора ":=" из языка «Паскаль».
Теперь можно построить гистограммы:
par(mfrow=c(1, 2))
hist(data, breaks=0:7, right=F, col="seagreen", main="Гистограмма НД1", xlab=x.txt, ylab="Число респондентов")
hist(data_ol, breaks=0:8, right=F, col="slateblue1", main="Гистограмма НД2", xlab=x.txt, ylab="Число респондентов")
Первая строка нужна, чтобы гистограммы вывелись рядом. Без нее вторая гистограмма затрет первую. Вот, что получилось:
Напоминает результат опроса, верно? Верно, особенность нашего исследования в том, что данные сгенерированы на основе гистограммы. Но данный шаг не лишен смысла, т.к.
- в ЖЖ нелинейный масштаб (скорее всего из-за количества голосов в первом варианте ответа);
- обе гистограммы изображены в одном масштабе и ориетированы вертикально, что позволяет уже сейчас производить сравнение с плотностью вероятности (probability density function) нормального распределения.
Следующий шаг имеет мало смысла для столь дискретизированного набора данных, как у нас. Он приведен здесь только для ознакомления с функцией density, которая строит более «сглаженную» (читай, усредненную) гистограмму по выборке.
plot(density(data), type="l", col="seagreen", lwd=2, main="Сглаженная плотность НД1")
plot(density(data_ol), type="l", col="slateblue1", lwd=2, main="Сглаженная плотность НД2")
Результат:
Вычислим выборочные параметры распределений.
mu = mean(data)
mu
mu_ol = mean(data_ol)
mu_ol
var(data)
var(data_ol)
sig = sd(data)
sig
sig_ol = sd(data_ol)
sig_ol
library(PerformanceAnalytics)
skewness(data)
skewness(data_ol)
kurtosis(data)# excess kurtosis (-3)
kurtosis(data_ol)
Результаты:
| № НД | Мат.ожидание | Дисперсия | Стандартное отклонение | Асимметрия (skewness) | Эксцесс (excess kurtosis) |
|---|---|---|---|---|---|
| 1 | 2.408437 | 1.708542 | 1.307112 | 0.4124443 | 0.1001578 |
| 2 | 2.465034 | 2.17858 | 1.476001 | 0.7198767 | 0.7943986 |
Как видно из таблицы, изменения во втором наборе данных
- едва ли изменили в среднем ожидаемое значение,
- увеличили разброс случайной величины,
- почти в 2 раза увеличили искаженность распределения вправо (у плотности распределения удлиннился правый «хвост»),
- почти в 8 раз увеличили толщину «хвостов» (по сравнению с нормальным распределением).
Сравним эмпирические функции распределения (ЭФР) с функциями распределения (cumulative distribution function) соответствующих нормальных распределений (N(2.408437, (1.307112)2) и N(2.465034, (1.476001)2)).
n1 = length(data)
plot(sort(data), (1:n1)/n1, type="S", col="seagreen", main="ЭФР НД1", xlab=x.txt, ylab="")
x = seq(0, 6, by=0.25)
lines(x, pnorm(x, mean=mu, sd=sig), type="l", col="orange", lwd=2)
n2 = length(data_ol)
plot(sort(data_ol), (1:n2)/n2, type="S", col="slateblue1", main="ЭФР НД2", xlab=x.txt, ylab="")
x2 = seq(0, 7, by=0.25)
lines(x2, pnorm(x2, mean=mu_ol, sd=sig_ol), type="l", col="orange", lwd=2)
Вывод:
От функций распределения перейдем к квантилям (quantile), обратным функциям распределения.
quantile(data)
quantile(data_ol)
qnorm(p=c(0, .25, .5, .75, 1), mean=mu, sd=sig)
qnorm(p=c(0, .25, .5, .75, 1), mean=mu_ol, sd=sig_ol)
В нашем конкретном случае этап довольно скучный, т.к. выборки отличаются только сотым процентилем:
| Распределение | q0 | q.25 | q.5 | q.75 | q1 |
|---|---|---|---|---|---|
| НД1 | 0 | 2 | 2 | 3 | 6 |
| N(2.408437, (1.307112)2) | -Inf | 1.526803 | 2.408437 | 3.290070 | Inf |
| НД2 | 0 | 2 | 2 | 3 | 7 |
| N(2.465034, (1.476001)2) | -Inf | 1.469486 | 2.465034 | 3.460582 | Inf |
И если НД1 квартилями походит на нормальное распределение хотя бы с округлением, то НД2 даже это не помогает.
Схема квантилей (normal Q-Q plot) для наших сильно дискретизированных выборок не сильно полезна. Упоминается, дабы осветить функцию qqnorm.
qqnorm((data-mu)/sig, col="seagreen")
abline( 0, 1, col="orange", lwd=2)
qqnorm((data_ol-mu_ol)/sig_ol, col="slateblue1")
abline( 0, 1, col="orange", lwd=2)
Результат выглядит не захватывающе, зато веселенько:
И завершает список наглядных выводов ящик с усами (boxplot).
boxplot(data, outchar=T, main="Ящик с усами НД1", ylab=x.txt)
boxplot(data_ol, outchar=T, main="Ящик с усами НД2", ylab=x.txt)
Графика:
Построение наглядно отражает робастные характеристики выборки (устойчивые к наличию выбросов):
- первый и третий квартиль (верхняя и нижняя граница прямоугольника),
- второй квартиль (медиана, утолщенная горизонтальная линией),
- доверительный интервал (верхний и нижний «ус», все за этими пределами считается выбросами),
- собственно, выбросы (окружности за пределами «усов»).
Доверительный интервал в данном случае считается примерно как отступ от первого/третьего квартиля на 1,5 интерквартильного размаха. За подробностями — ?boxplot.
Заключение
НД1 отклоняется меньше НД2 от нормального распределения в виду:
- меньших значений асимметрии и эксцесса,
- меньшего различия квантилей распределения с выборочным ожиданием и выборочным стандартным отклонением в сравнении с квантилями выборки.
Дополнительная информация
Альтернативные вводные материалы в R (англ.):
Вторая и третья ссылки — часть официальной документации. Если есть ссылки на дельные вводные статьи на великом и могучем, пишите — добавлю.
Основной целью статьи является привлечение внимания общественности к R как инструменту анализа. Если кто-либо из знающих людей представит более углубленный материал, буду искренне рад и с удовольствием ознакомлюсь.
Корректоры — в личку. Остальные — добро пожаловать в комментарии.
Спасибо всем за внимание.
Автор: theoden
