Введение в R-project

^[1]Во всем Хабре сыскалась лишь пара ^[2] статей ^[3] на вышеуказанную тему. А тема благодатная. Да и в минувшую среду как раз окончился курс "Introduction to Computational Finance and Financial Econometrics ^[4]". По мотивам его пятой недели «Descriptive statistics» и появился этот пост. Причастившимся будет неинтересно, а желающих познакомиться с базовыми приемами анализа данных при помощи R  — прошу под хабракат.

Предварительные соглашения

О терминах

У автора из статистики был только семестр «тервера» N лет назад. Поэтому после сомнительно переведенных слов и их сочетаний будет указан исходный английский термин (курсивом в скобках). Специалисты, пожалуйста, шлите в личку более корректные варианты терминов. Спасибо.

Об установке

На установке ПО внимание не заостряется намеренно, в виду тривиальности. По крайней мере на Windows платформе все свелось к стандартному «далее -> далее ->… -> готово». Единственный требуемый для выполняемого в статье кода пакет PerformanceAnalytics устанавливается через меню «Пакеты / Установить пакет(ы)...», выбор ближайшего к вам зеркала, выбор нужного пакета из списка.

Набор данных

Хотелось избежать типичности: продажи, квартиры, рентабельность акций (simple returns), — сколько можно? Поэтому предметная область нашей выборки — вечна ^[5] как в контексте Хабра, так и вне его контекста. Не так давно в блоге СамиЗнаетеКого был опубликован опрос «Какого размера у вас грудь?» ^[6]. Учитывая, что в него были включены два варианта ответа для отсева нерелевантной аудитории, есть некоторая уверенность в правдоподобии выборки. Для удобства результаты приведены и здесь:

Цель

В рамках нашего мини-исследования сравним с нормальным распределением 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) вернет вектор из ста восьмидесяти четырех нулей. В рекомендациях гугла ^[7] и еще в нескольких источниках встречалось мнение, что символ равенства негоже использовать для присваивания, лучше — "<-". Знающие люди, пожалуйста, изложите в комментариях достаточно веские обоснования, чтобы писать 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="Число респондентов")

Первая строка нужна, чтобы гистограммы вывелись рядом. Без нее вторая гистограмма затрет первую. Вот, что получилось:

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

1. в ЖЖ нелинейный масштаб (скорее всего из-за количества голосов в первом варианте ответа);
2. обе гистограммы изображены в одном масштабе и ориетированы вертикально, что позволяет уже сейчас производить сравнение с плотностью вероятности (probability density function) нормального распределения ^[8].

Следующий шаг имеет мало смысла для столь дискретизированного набора данных, как у нас. Он приведен здесь только для ознакомления с функцией 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)

Результаты:

Как видно из таблицы, изменения во втором наборе данных

• едва ли изменили в среднем ожидаемое значение,
• увеличили разброс случайной величины,
• почти в 2 раза увеличили искаженность распределения вправо (у плотности распределения удлиннился правый «хвост»),
• почти в 8 раз увеличили толщину «хвостов» (по сравнению с нормальным распределением).

Сравним эмпирические функции распределения (ЭФР) с функциями распределения (cumulative distribution function) соответствующих нормальных распределений (N(2.408437, (1.307112)²) и N(2.465034, (1.476001)²)).

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)

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

И если НД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)

Результат выглядит не захватывающе, зато веселенько:

И завершает список наглядных выводов ящик с усами ^[9] (boxplot).

boxplot(data, outchar=T, main="Ящик с усами НД1", ylab=x.txt)
boxplot(data_ol, outchar=T, main="Ящик с усами НД2", ylab=x.txt)

Графика:

Построение наглядно отражает робастные ^[10] характеристики выборки (устойчивые к наличию выбросов ^[11]):

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

Доверительный интервал в данном случае считается примерно как отступ от первого/третьего квартиля на 1,5 интерквартильного размаха. За подробностями — ?boxplot.

Заключение

НД1 отклоняется меньше НД2 от нормального распределения в виду:

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

Дополнительная информация

Альтернативные вводные материалы в R (англ.):

• R Introduction ^[12], сопроводительные скрипты ^[13]
• An Introduction to R ^[14]
• R for Beginners ^[15]

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

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

Корректоры — в личку. Остальные — добро пожаловать в комментарии.
Спасибо всем за внимание.

Автор: theoden

Источник ^[16]