Статистичне програмування на R: Частина 2. Функціональне програмування та аналіз даних (исходники), Різне, Програмування, статті

R – потужний функціональний мова програмування і середу аналізу наборів статистичних даних. Будучи середовищем аналізу, R дозволяє створювати різні графічні представлення даних з командного рядка. Читайте першу частину цієї серії статей і статтю Камерона Лейрда про R – з них можна більше дізнатися про R, про платформах, на яких може працювати R, а також про набір доступних для неї пакетів.


Як і інтерактивні оболонки Python, Ruby, wish для Tcl / Tk та багато середовища Lisp, оболонка R відмінно підходить для дослідження операцій – використання редагованих команд дозволяє здійснювати різну обробку даних. На відміну від інтерактивних оболонок багатьох інших мов програмування, але в повній відповідності з орієнтацією R на дані, оболонка R дозволяє зберігати всю середу проекту (по одній на кожен робочий каталог). Одна з корисних особливостей – Історія команд: її можна переглядати, редагувати і зберігати у файлі. Rhistory. Проте основна особливість збереження середовища робочих проектів полягає в тому, що всі дані зберігаються в бінарній формі у файлі. RData. Радимо регулярно використовувати команду rm(), Щоб зберігається середу даних не розросталася необмежено (remove() – Синонім для rm()).


Ця стаття – друга частина розповіді про R. Перша стаття знайомить читача з деякими базовими типами даних R, включаючи вектори і багатовимірні масиви (2-мірні масиви ще називають матрицями), масиви даних для “розумних” таблиць, списки для гетерогенних колекцій і т.д. У минулій статті були також виконані базовий статистичний аналіз і побудова графіків для великого масиву даних, зібраного одним із співавторів (Бредом): ці дані – річна історія зміни температури навколо будинку Бреда, зібрана з трихвилинним інтервалами часу. Як і більшість реальних даних, ці температурні дані містять дефекти і похибки виміру – і в цій статті ми будемо намагатися виявити підозрілі дані.


Наша стаття вирішує два завдання. По-перше, ми продовжимо аналіз даних, які розглядалися в попередній статті, щоб детальніше розглянути можливості самої мови R. По-друге, ми досліджує загальні закономірності даних і покажемо, як знаходити їх, використовуючи R.


R як мова програмування


Розповсюджуваний за ліцензією GPL мову програмування R має двох предків: S / S-PLUS і Scheme (Lisp). Особливо слід підкреслити зв’язок зі Scheme, так як основні особливості функціонального програмування в R успадковані від Scheme. Уважні читачі могли відзначити визначна відсутність в попередній статті одного істотного аспекту – управління виконанням програми! У цій статті ми як і раніше не зачіпаємо цю тему. Насправді в R є відмінні команди if, else, while і for, Дуже схожі на аналогічні команди Python (в інших мовах синтаксис цих команд дещо інший). Плюс до того є ще й команди repeat, break, next і switch. Програмісти, які використовують процедурні мови, будуть здивовані тим, наскільки багато можна зробити без управління виконанням – і ще більше вони будуть збентежені, дізнавшись, наскільки простіше стають без цього багато завдань!


Оголошуйте те, що вам потрібно


Цілком у дусі функціонального програмування в R можна зробити практично все за допомогою простих декларативних операторів. В R є дві особливості, які в більшості випадків роблять управління виконанням зайвим. По-перше, ви вже бачили, що більшість операцій над об’єктами колекцій працюють поелементно. Щоб зробити щось над елементами вектора даних, не потрібно цикл, тому що можна просто виконати операцію відразу з усіма елементами вектора:


Лістинг 1. Поелементні операції в R




> A = 1:10 # Діапазон значень від 1 до 10 > B = a +5 # Додаток 5 до кожного елементу > B # Висновок вектора b
[1] 6 7 8 9 10 11 12 13 14 15

Також можна оперувати тільки елементами з певними індексами, використовуючи “індексний масив”:


Лістинг 2. Використання індексних масивів для вибору елементів




> C = b # Копіюємо b в a > Pos = c (1,2,3,5,7) # Змінюємо первинні індекси > C [pos] = c [pos] * 10 # перерозподіляємо позиції індексів > C # Висновок вектора
[1] 60 70 80 9 100 11 120 13 14 15

Або, що, можливо, краще за все, можна використовувати синтаксис, схожий на спискові висловлювання мов Haskell або Python, і оперувати тільки над елементами, які мають необхідну властивість:


Лістинг 3. Використання предикативного вибору елементів




> d = c > D [d%% 2 == 0] = -1 # Присвоєння -1 всім парних елементів
> d
[1] -1 -1 -1 9 -1 11 -1 13 -1 15

Уважні читачі можуть зауважити, що в цих прикладах використовується присвоювання за допомогою знака =, Тоді як в більшості прикладів раніше використовувалися знаки <-. Знак рівності робить те ж, що й стрілка, але може бути використаний тільки в високорівневих операціях, але не в у вкладених виразах. Якщо є сумніви з приводу використання знака рівності, безпечніше використовувати стрілку.


Аналізуємо дані про температуру


У першій статті ми за допомогою функції read.table() вважали приблизно 170 000 свідчень кожного з 4 датчиків. Однак для більш зручного доступу до окремих серій вимірювань ми скопіювали колонки даних в вектори, названі за місцями розміщення датчиків: outside (на вулиці), livingroom (вітальня), basement (підвал) і lab (лабораторія). Не потрібно забувати, що деякі результати вимірювань були загублені: іноді на записувальному комп’ютері відбувалися збої, прилади іноді відмовляли, не всі 4 термометра були включені точно в один день і тощо. Іншими словами, наші температурні дані багато в чому схожі на реальні дані, з якими доводиться мати справу в житті.


При первинному дослідженні температурних даних була помічена аномалія в гистограммах. Вона проявляється великим піком зовнішньої температури поблизу 24 градусів Цельсія. Наше перше припущення полягало в тому, що ці піки відображають розсинхронізація при записи, коли температура з внутрішнього приміщення записувалася в дані зовнішнього приміщення (і, відповідно, навпаки). Щоб проілюструвати можливості R, подивимося, чи можна довести або спростувати це припущення.


Пошук аномалій


Перший крок у дослідженні аномалій – це їх знаходження. А саме, точкова аномалія повинна характеризуватися великими стрибками температури з обох сторін від точки даних. Ще конкретніше, ми можемо очікувати, що 2 сусідні з аномалією точки будуть звичайно набагато вище або набагато нижче, ніж точка потенційної аномалії. Наприклад, проста послідовність температур (з трихвилинним інтервалом) “20, 16, 13” демонструє незвично швидке падіння температури, але не породжує підозри на одноточечную помилку в середньому відліку. Звичайно, не можна апріорі виключити існування інших типів помилок або похибок, які зачіпають не тільки одиничні відліки даних.


Наші перші ідеї щодо ідентифікації неправильних вимірювань були досить складні. Можна звернути увагу на високочастотні компоненти перетворення Фур’є наших послідовностей. Або можна взяти похідні векторів температури (можливо, другі похідні, щоб знайти вигини графіка). Можна спробувати виконати згортки температурних векторів. Всі ці операції вбудовані в R. Але потім ми залишили високі помисли і спустилися на землю. Ми шукаємо просту закономірність – це відліки температур, що мають великі скалярні різниці з сусідніми вимірами. Іншими словами, нам достатньо простого віднімання.


Щоб знайти всі різниці між сусідніми точками даних, ми створимо таблицю даних, колонки якої будуть відповідати попереднім, поточним і наступним точкам даних. Це дозволить нам фільтрувати всю таблицю даних, вибираючи цікавлять ряди.


Лістинг 4. Пошук одиночних похибок в даних вуличної температури




> Len <- length (outside) # Короткий ім'я довжини вектора > I <- 1: (len-2) # Короткий ім'я вектора номерів рядків колонок таблиці даних > # Створення таблиці даних для вікна вимірювань на кожен рядок
> odf <- data.frame(lst=outside[1:(len-2)],
+ now=outside[2:(len-1)],
+ nxt=outside[3:len] ) > # Створення вектора локальних змін температури, додавання до таблиці даних
> odf$flux <- (odf[i,”now”]*2) – (odf[i,”lst”] + odf[i,”nxt”]) > Odf2 <- odf [! Is.na (odf $ flux),] # Виключення "Not Available" вимірів > Oddities <- odf2 [abs (odf2 $ flux)> 6,] # Дивно, якщо відхилення більше 6

Отже, що ми маємо після фільтрації? Давайте подивимося:


Лістинг 5. Перегляд похибок в даних вуличної температури




> oddities
lst now nxt flux
2866 21.3 15.0 14.9 -6.2
79501 -1.5 -6.2 -4.1 -6.8
117050 21.2 24.6 21.6 6.4
117059 20.6 23.4 20.1 6.1
127669 24.1 21.2 24.7 -6.4
127670 21.2 24.7 21.5 6.7

Ці дані говорять про те, що є точки з сильним відзнакою щодо сусідніх вимірів. Але більша частина з них не виглядає як кандидати на помилки синхронізації. Наприклад, в інтервалі часу 79501 температура 6.2 оточена двома явно вищими температурами. Однак жодна з них не виглядає схожою на температуру всередині будинку, швидше вже це дуже холодний порив вітру.


З іншого боку, схоже, що ми маємо перестановки близько тимчасових інтервалів 117059 і 12670. Середні значення температур близькі до внутрішньої (24 градуса), а сусідні вимірювання значно нижче – хоча це і не мороз. Може бути, це вимірювання, зроблені навесні.


Оформлення корисних операцій у функцію


Тепер нам потрібно дізнатися, чи мають наші кандидати на перестановку відповідності в інших відліку. Якщо ми знайдемо відповідні дані вуличної температури в одному з наборів даних температур всередині будинку, це підтвердить нашу гіпотезу. Але ми не хочемо заново набирати програму, змінюючи скрізь outside на lab. Що дійсно варто зробити – це об’єднати всю послідовність операцій в одну функцію:


Лістинг 6. Оформлення процесу виявлення відхилень у функцію




oddities <- function(temps, flux) {
len <- length(temps)
i <- 1:(len-2)
df <- data.frame(lst=temps[1:(len-2)],
now=temps[2:(len-1)],
nxt=temps[3:len])
df$flux <- (df[i,”now”]*2) – (df[i,”lst”] + df[i,”nxt”])
df2 <- df[!is.na(df$flux),]
oddities <- df2[abs(df2$flux) > flux,]
return(oddities)
}

Маючи функцію, значно простіше обробляти різні набори даних і встановлювати пороги відхилень. Запуск функції oddities(lab,6)не видає відміток часу наших кандидатів на перестановку (так само як і для livingroom і basement). Проте погляд на записи температур лабораторії підносить інший сюрприз:


Лістинг 7. Величезні температурні коливання в лабораторії




> oddities(lab, 30)
lst now nxt flux
47063 19.9 -2.6 19.5 -44.6
47268 17.7 -2.6 18.2 -41.1
84847 17.1 -0.1 17.0 -34.3
86282 14.9 -1.0 14.8 -31.7
93307 14.2 -6.4 14.1 -41.1

Показань такого сорту ми не очікували. Може бути, краще пояснення – це те, що Бред відкривав вікно лабораторії в жорстокий мороз. Якщо так, то Девід уражається ефективності каміна Бреда, відновлюючого тепло всього за 3 хвилини вентиляції кімнати.


Правильне пояснення


Можливо, перегляд повних даних по цікавлять нас моментів часу роз’яснить ситуацію:


Лістинг 8. Повний набір записів температур поблизу кроку 47063




> glarp[47062:47066,]
timestamp basement lab livingroom outside
47062 2003-10-31T17:07 21.5 20.3 21.8 -2.8
47063 2003-10-31T17:10 21.3 19.9 21.2 -2.7
47064 2003-10-31T17:13 20.9 -2.6 20.9 -2.6
47065 2003-10-31T17:16 20.8 19.5 20.8 -2.6
47066 2003-10-31T17:19 20.5 19.4 20.7 -2.8

По-перше, зауважимо, що позначки часу, які повертаються функцією oddities(), Завищені на одиницю. Ах так, ми використовували зсув для створення вікна для кожного рядка таблиці даних. При цьому, вся сукупність даних свідчить на користь ідеї “відкритих дверей” – напередодні Хеллоуїна в 5:00 вечора в Колорадо може бути досить прохолодно, а температура може сильно змінюватися, коли в двері Бреда заглядають гуляють (і отримують свої цукерки за 3 хвилини). Отже, можливо, ця загадка дозволена.


І все-таки, що ж щодо загадки, що почала наші дослідження? Чому з’являються піки близько 24 градусів у вимірі вуличної температури? Може бути, необхідно подивитися на гістограму більш уважно? Зокрема, можна використовувати предикативне критерій для індексування вектора і звузити гістограму до даного нас температурного діапазону:


Лістинг 9. Звуження температурного діапазону




hist(outside[outside < 26 & outside > 23],
breaks=90, col=”green” border=”blue”)

Рисунок 1. Гістограма для вулиці в звуженої температурної області
Гістограма для вулиці в звуженої температурної області

Подивившись докладніше на локалізовані температурні піки, можна виявити окремі западини праворуч від кожного піку. Вони видні, якщо подивитися на відхилення в десяті частки для проміжку поруч зі значенням 24.7 градуса. Найбільш ймовірно, що це є наслідком похибки округлення термометра. Відмінно, але все ж ми ще не впевнені. Навіть після деякого згладжування невеликий пік таки залишається.


Проміжний статистичний аналіз


Одна із сильних сторін R – це здатність розраховувати лінійні і нелінійні моделі регресій. Розглянемо простий приклад. Створимо два вектори: x буде часом в днях, починаючи з дати початку збору даних, а y буде відповідною температурою на вулиці.


Лістинг 10. Створення вектора регресії




> y <- glarp$outside
> x <- 1:length(y)/(24*60/3)

Можна підібрати пряму для цих значень:


Лістинг 11. Підбір лінійної регресії




> l1 = lm(y ~ x)
> summary(l1)
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-29.4402 -7.4330 0.2871 7.4971 23.1355
Coefficients:
Estimate Std. Error t value Pr(>/t/)
(Intercept) 10.2511014 0.0447748 228.95 <2e-16 ***
x -0.0037324 0.0002172 -17.19 <2e-16 ***

Signif. codes: 0 `***” 0.001 `**” 0.01 `*” 0.05 `.” 0.1 ` ” 1
Residual standard error: 9.236 on 169489 degrees of freedom
Multiple R-Squared: 0.00174, Adjusted R-squared: 0.001734
F-statistic: 295.4 on 1 and 169489 DF, p-value: < 2.2e-16

Знак “~” означає формулу. Ця дія вказує R знайти коефіцієнти A і B, які мінімізують формулу sum((y[i]-(A*x[i]+B))^2). Кращі значення – коли A одно -0.0037324 (дуже близький до нуля кут), а B одно 10.2511014. Зауважимо, що залишкова стандартна похибка дорівнює 9.236, що мало відрізняється від стандартного відхилення самого y. Це говорить про те, що проста лінійна функція від часу – дуже погана модель для аналізу вуличної температури.


Можливо, краще було б використати апроксимацію синусами і косинусами з періодами 1 день і 1 рік. Таку модель можна спробувати, змінивши формулу на наступну:


Лістинг 12. Спроба наближення тригонометричними кривими




> l2 = lm(y ~ +I(sin(2*pi*x/365)) +I(cos(2*pi*x/365))
+ +I(sin(2*pi*x)) +I(cos(2*pi*x)) )

Синтаксис цієї формули дуже хитрий: всередині викликів I() арифметичні функції мають своє звичайне призначення. Наприклад, у перший виклик I() записано такий вираз: 2 помножити на pi, помножити на x і ділити на 365. Знак “+” перед I() вказує не на додавання, а на те, що вираз після + повинно бути включено в модель. Результат може бути відображений, як звичайно, за допомогою команди summary().


Лістинг 13. Результати тригонометричної регресії




> summary(l2)
Call:
lm(formula = y ~ +I(sin(2 * pi * x/365)) +I(cos(2 * pi * x/365))
+I(sin(2 * pi * x)) +I(cos(2 * pi * x)))
Residuals:
Min 1Q Median 3Q Max
-21.7522 -3.4440 0.1651 3.7004 17.0517
Coefficients:
Estimate Std. Error t value Pr(>/t/)
(Intercept) 9.76817 0.01306 747.66 <2e-16 ***
I(sin(2 * pi * x/365)) -1.17171 0.01827 -64.13 <2e-16 ***
I(cos(2 * pi * x/365)) 10.04347 0.01869 537.46 <2e-16 ***
I(sin(2 * pi * x)) -0.58321 0.01846 -31.59 <2e-16 ***
I(cos(2 * pi * x)) 3.64653 0.01848 197.30 <2e-16 ***

Signif. codes: 0 `***” 0.001 `**” 0.01 `*” 0.05 `.” 0.1 ` ” 1
Residual standard error: 5.377 on 169486 degrees of freedom
Multiple R-Squared: 0.6617, Adjusted R-squared: 0.6617
F-statistic: 8.286e+04 on 4 and 169486 DF, p-value: < 2.2e-16

Залишкова помилка все ще велика (5.377), але утішимо себе тим фактом, що погода в Колорадо виключно непередбачувана.


R надає безліч інструментів для аналізу тимчасових даних. Наприклад, ми можемо побудувати графік функції автокореляції для температури у вітальні:


Лістинг 14. Функція автокореляції температури у вітальні




> acf(ts(glarp$livingroom, frequency=(24*60/3)),
+ na.action=na.pass, lag.max=9*(24*60/3))

Рисунок 2. Графік функції автокорелляціі температури у вітальні

Вбудована функція ts() створює тимчасової ряд з вектора glarp$livingroom. Частота вибірки задається в відліку в день. Не дивно, що температура сильно корелює, коли зсув становить цілу кількість днів. Помітні також невеликі піки на семи днях. Причина цього в те, що термостат Бреда встановлює іншу температуру на вихідних (коли Бред зазвичай йде з дому на весь день), в результаті, кореляції при установці семиденного вікна виявляються трохи більше.


Висновок


Отже, ми використовували R для аналізу структури та аномалій в наборах даних, які мають точно такі ж недоліки та потенційні проблеми, як і практично всі реальні наукові дані. У цій статті ми також побачили, як стиль функціонального програмування, реалізований в R, допомагає швидко досліджувати закономірності, дані та аналітичні сценарії. У третій частині цієї серії ми продовжимо пошук закономірностей у великих наборах даних, використовуючи більш складні статистичні методи (але все одно зачіпаючи лише малу частину найбагатших можливостей R).


Схожі статті:


Сподобалася стаття? Ви можете залишити відгук або підписатися на RSS , щоб автоматично отримувати інформацію про нові статтях.

Коментарів поки що немає.

Ваш отзыв

Поділ на параграфи відбувається автоматично, адреса електронної пошти ніколи не буде опублікований, допустимий HTML: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>

*

*