Статистическое программирование на R: Часть 1. Купаемся в изобилии статистических возможностей

Приведение исходных данных к единому виду

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

Девид Мертц (David Mertz), Разработчик, Gnosis Software, Inc.

Девид Мертц считает, что искусственные языки вполне естественны, а естественные кажутся немного искусственными. Вы можете связаться с ним по mertz@gnosis.cx; вы можете исследовать все стороны его жизни на его личной web-странице. Посмотрите его книгу, Text Processing in Python. Пожелания и предложения по поводу прошлых и будущих статей только приветствуются.



Бред Хантинг, Mathematician, Университет Колорадо

Бред занимается сетями и системным администрированием UNIX®в различных компаниях уже 14 лет. В настоящий момент он планирует стать кандидатом наук по прикладной математике в университете Колорадо в Болдере, и оплачивает свои счета, осуществляя поддержку UNIX в отделе Computer Science.



24.07.2008

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

Разработчики описывают цели создания R в статье "An Introduction to R" ("Введение в R"; см. ссылку на полный текст статьи ниже в разделе Ресурсы):

Рекомендуется использовать для различных проектов в R отдельные рабочие каталоги. Достаточно часто во время работы создаются объекты с именами типа x или y. Такие имена имеют смысл в контексте определенного проекта, однако будет сложно определить, что имелось в виду, когда [...] файлы нескольких проектов сохранялись в одном каталоге.

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

Распространяемый по лицензии GPL язык программирования R имеет двух родителей: коммерческий язык программирования S/S-PLUS, из которого была позаимствована большая часть синтаксиса, и язык программирования Scheme, предоставивший многие (наиболее изящные) семантические структуры. Ранние реализации S датируются 1984 годом; более поздние версии (включая S-PLUS) добавили много новых возможностей. Scheme (как и Lisp), конечно, появился еще раньше. R возник как свободная реализация расширенной версии S/S-frPLUS в 1997 году и сразу стал пользоваться успехом у пользователей и разработчиков. Не сомневайтесь в полезности R. Камерон Лейрд представил ряд материалов и ссылок, касающихся R, в статье на developerWorks R handy for crunching data.

Среда R доступна в бинарном виде для многих компьютерных платформ: Linux, Windows, Mac OS X и Mac OS Classic. Естественно, предоставляется также исходный код для компиляции на других платформах (например, соавтор этой статьи Брэд без труда скомпилировал R на FreeBSD). Работа R не всегда корректна на некоторых платформах: например, окна графиков, использующих Quartz, на компьютере Дэвида под Mac OS X зависают; более того, на компьютере Брэда (FreeBSD/AMD Athlon) выход из R может спровоцировать перезагрузку (возможно, это происходит из-за некорректных опций ядра SSE, но такое поведение все равно неприятно). Несмотря на это, R в целом стабилен, быстр и обладает совершенно удивительным набором статистических и математических функций. Вдобавок к огромному набору стандартных пакетов и функций существуют дополнительные пакеты.

Модель данных R

Базовым объектом данных в R является вектор. Различные виды векторов добавляют такие возможности, как работа с (многомерными) массивами, структурами данных, (разнородными) списками и матрицами. Так же как и в NumPy/NumArray или Matlab, операции над векторами и однородными элементами производятся поэлементно. Несколько простых примеров работы в R дадут возможность понять его синтаксис (в примеры включены команды и ответы системы):

Листинг 1. Векторы и элементарные операции
> a <- c(3.1, 4.2, 2.7, 4.1)  # Присваивание «стрелкой» предпочтительнее, чем "="
> c(3.3, 3.4, 3.8) -> b       # Можно также присваивать направо
> assign("c", c(a, 4.0, b))      # или явно назначать имя переменной
> c                              # Конкатенация "уплощает" аргументы
[1] 3.1 4.2 2.7 4.1 4.0 3.3 3.4 3.8
> 1/c                            # Операция над каждым элементом вектора
[1] 0.3225806 0.2380952 0.3703704 0.2439024 0.2500000 0.3030303 0.2941176
[8] 0.2631579
> a * b                          # Цикл по более короткому вектору "b" (с предупреждением)
[1] 10.23 14.28 10.26 13.53
Warning message:
longer object length
        is not a multiple of shorter object length in: a * b
> a+1                            # "1" обрабатывается как вектор длины 1
[1] 4.1 5.2 3.7 5.1

Ниже приведены другие примеры индексации, разделения, именованных и необязательных аргументов и других элементов синтаксиса R. Командная оболочка R, особенно если у вас установлен GNU readline, -- это замечательный интерфейс для обучения. Заведите привычку пользоваться командой help(function), чтобы учиться в процессе работы (также можно использовать команду ?function). Пользователи оболочки Python найдут в оболочке R много знакомого и полезного.


Набор данных о температуре

Брэд почти в течение года собирал данные о температуре с четырех термометров внутри и вокруг своего дома и с помощью GnuPlot автоматически преобразовывал изменения температуры в скользящие по времени графики, доступные через Web (более подробную информацию о GnuPlot см. в разделе Ресурсы).Хотя такая любительская коллекция данных не имеет реальной научной ценности, она имеет ряд замечательных особенностей, которые делают ее похожей на настоящие научные данные.

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

Чтение данных

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

Листинг 2. Формат файла исходных данных о температуре
2003 07 25 16 04 27.500000
2003 07 25 16 07 27.300000
2003 07 25 16 10 27.300000

Изначально данные могут быть прочитаны, например, следующим образом:

Листинг 3. Первоначальное чтение данных о температуре
> lab <- read.fwf('therm/lab', width=c(17,9))   # Формат фиксированной ширины
> basement <- read.fwf('therm/basement', width=c(17,9))
> livingroom <- read.fwf('therm/livingroom', width=c(17,9))
> outside <- read.fwf('therm/outside', width=c(17,9))
> l_range <- range(lab[,2]) 	  # Вектор блока данных: весь второй столбец
> b_range <- range(basement[,2])    # range() дает min/max
> v_range <- range(livingroom[,2])
> o_range <- range(outside[,2])
> global <- range(b_range, l_range, v_range, o_range)
> global          			  # Интервалы температур по всем точкам
[1] -19.8  32.2

Очистка данных

Исходный формат данных имеет некоторые проблемы. Например, не указываются потерянные данные - вместо этого стоят простые метки отсутствия строчек и отметки времени. Кроме того, даты хранятся в нестандартном формате (отличном от ISO8601/W3C). Менее существенный недостаток - избыточное повторение временных меток в каждом из четырех файлов. Конечно, мы можем очистить данные в самом R, но вместо этого воспользуемся рекомендациями авторов R, представленными в документе "R Data Import/Export" («Импорт/Экспорт данных в R», см. раздел Ресурсы). Обработку текстовых данных обычно лучше производить на языке, предназначенном для этого: в доказательство мы написали скрипт на языке Python, создающий единый файл данных, который будет напрямую загружаться в R. Для примера приведем несколько первых строк нового файла данных glarp.temps:

Листинг 4. Формат объединенных данных о температуре
Время            Подвал    Лаборатория Гостиная   Улица
2003-07-25T16:04 24.000000      NA     29.800000  27.500000
2003-07-25T16:07 24.000000      NA     29.800000  27.300000
2003-07-25T16:10 24.000000      NA     29.800000  27.300000

Теперь будем работать с обновленными данными:

Листинг 5. Работа с объединенными температурными данными
> glarp <- read.table('glarp.temps', header=TRUE, as.is=TRUE)
> timestamps <- strptime(glarp[,1], format="%Y-%m-%dT%H:%M")
> names(glarp)              # Каковы имена колонок?
[1] "timestamp"  "basement"   "lab"        "livingroom" "outside"
> class(glarp[,'basement'])  # Какой вид данных в колонке «подвал»?
[1] "numeric"
> basement <- glarp[,2]         # присваивание по позиции
> lab <- glarp[,'lab']          # присваивание по имени
> outside <- glarp$outside      # соответствует предыдущему присваиванию
> livingroom <- glarp$living    # уникальное исходное имя
> summary(glarp)  		
# Удобное встроенное средство описания большинства объектов в R
  timestamp            basement            lab            livingroom
 Length:171349      Min.   :   6.40   Min.   :  -6.40   Min.   :   7.20
 Class :character   1st Qu.:  17.00   1st Qu.:  16.60   1st Qu.:  18.10
 Mode  :character   Median :  19.10   Median :  17.90   Median :  20.30
                    Mean   :  18.88   Mean   :  18.12   Mean   :  20.17
                    3rd Qu.:  20.50   3rd Qu.:  19.50   3rd Qu.:  22.00
                    Max.   :  27.50   Max.   :  25.50   Max.   :  31.30
                    NA's   :1854.00   NA's   :2406.00   NA's   :1855.00
    outside
 Min.   : -19.800
 1st Qu.:   2.100
 Median :   9.800
 Mean   :   9.585
 3rd Qu.:  17.000
 Max.   :  32.200
 NA's   :1858.000

Базовый статистический анализ

Взгляните на функцию range(): min() и max() находят отдельные экстремумы интервала данных. summary(), очевидно, показывает ту же информацию, но в виде, недоступном для прямого использования в других вычислениях. Начнем с нахождения нескольких базовых статистических свойств этих данных:

Листинг 6. Базовые статистические расчеты
> mean(basement)        	# Среднее значение не будет вычислено, если включить
[1] NA
> mean(basement, na.rm=TRUE)
[1] 18.87542
> sd(basement, na.rm=TRUE)  	
# При вычислении среднего отклонения также нужно исключить отсутствующие данные
[1] 2.472855
> cor(basement, livingroom, use="all.obs")   # Все наблюдения: не проходит
Error in cor(basement, livingroom, use = "all.obs") :
        missing observations in cov/cor
> cor(basement, livingroom, use="complete.obs")
[1] 0.9513366
> cor(outside, livingroom, use="complete.obs")
[1] 0.6446673

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

Распределение температур

Мы посчитали среднее значение и среднее отклонение температур; интуитивно можно полагать, что температуры будут распределены по нормальному закону. Давайте проверим:

Листинг 7. Создание гистограммы одной короткой строкой
> hist(outside)

Многие команды R приводят к появлению всплывающих окон с графиками, чертежами или диаграммами данных. Детали создания всплывающих окон зависят от платформы и настроек. Графики можно также перенаправить для сохранения во внешних файлах для дальнейшего использования. Команда hist() приводит к следующему результату:

Рисунок 1.
Гистограмма температур на улице

Неплохо для первой попытки. С помощью нескольких параметров можно сузить интервал усреднения:

Листинг 8. Сужение диапазона суммирования гистограммы
hist(outside, breaks=1000, border="blue")
Рисунок 2.
Ужатая гистограмма температур на улице

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

Еще пара небольших изменений и распределение внутренней температуры:

Листинг 9. Гистограммы распределения температур в гостиной
> hist(livingroom, breaks=40, col="blue", border="red")
> hist(livingroom, breaks=400, border="red")
Рисунок 3.
Гистограмма температур в гостиной с шагом 40
Рисунок 4.
Гистограмма температур в гостиной с шагом 400

График распределения температур в гостиной представляется более разумным. Сильные неравномерности, проявляющиеся при высоком разрешении, похоже, являются следствием небольших инструментальных погрешностей. Но основное трехпиковое распределение -- это результат работы управляемого таймером термостата Брэда (большой пик в районе 21, меньшие -- в районе 16 и 24 градусов).


Подробнее о визуализации данных

Каждому пункту измерений соответствует линейный вектор значений температур. Можно ожидать наличия в данных двух основных циклов: суточного и годичного (ночи и зимы холоднее).

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

Листинг 10. Преобразование вектора и создание графика температуры
> oarray <- outside[1:170880] # Необходимо отбросить несколько точек данных, 
относящихся к последнему дню
> dim(oarray) <- c(480,356)   # Преобразование вектора в двумерный массив
> plot(oarray[1,], col="blue", type="l", main="4 p.m. outside temp",
+      xlab="Day of year (starting July 25, 2003)",
+      ylab="Temperature (Celsius)")
Рисунок 5.
Температура на улице в 4 вечера

Когда мы преобразовали данные в матрицу "время-день" (двумерный массив), естественный шаг -- взять значения температуры для каждого дня и построить график годичного поведения. Можно сделать иначе -- извлечь каждую 480-ую точку из вектора, но способ R более элегантен.

Трехмерный массив данных

Как насчет представления результатов измерения температуры за весь год? Один из подходов — цветокодированный график температур:

Листинг 11. Создание графика температур
> x <- 1:480                  # Создание индексов оси X
> y <- 1:356                  # Создание индексов оси Y
> z <- oarray[x,y]            # Определяем ось z (в действительности совпадающую с oarray)
> mycolors <- c(heat.colors(33),topo.colors(21))[54:1]
> image(x,y,z, col=mycolors,
+       main="Outside temperatures near Glarp", # Brad's name for his house
+       xlab="Minutes past 4 p.m",
+       ylab="Days past July 25, 2003" )
> dev2bitmap(file="outside-topo.pdf", type="pdfwrite")
Рисунок 6.
Карта температур в день х за весь год

Несколько комментариев о том, что было сделано. Определение осей и точек очевидно, если вы узнали нотацию создания списка чисел, подобную используемой в языке Python. Индексирование по x и y при создании z приводит к появлению массива с длиной и шириной индексов. В данном случае z тривиальна -- совпадает с oarray; однако достаточно легко систематично изменять значения или смещения, используемые для их получения. Цветовая карта mycolor подбирается путем проб и ошибок: мы полагаем, что необходимо использовать красные и желтые цвета для «горячих» температур (другими словами -- выше нуля градусов Цельсия) и неправильно использовать их для «холодных» температур. Поэтому мы добавляем в цветовой вектор несколько голубых/зеленых цветов. Мы также хотим использовать обратный порядок цветов вместо предлагаемого стандартными функциями цветовых карт -- это достаточно просто при использовании индексирования.

Вы можете заметить, что карта температур нарисована немного более контрастно, чем прежние графики. Команда image() создает на экране похожую, но менее впечатляющую картинку. Экспорт «текущего изображения» во внешний файл часто дает лучший результат, что и показано здесь.

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

Листинг 12. Создание графика топографической поверхности
> persp(x,y,z, theta=10, phi=60, ltheta=40, lphi=30, shade=.1, border=NA,
+       col=mycolors[round(z+20)], d=.5,
+       main="Outside temperatures near Glarp",
+       xlab="Minutes past 4 p.m",
+       ylab="Days past July 25, 2003",
+       zlab="Temperature", )
Рисунок 7.
Топографическая карта температур в день х за весь год

Заключение

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

Мы надеемся, что эта статья дала читателю представление о работе в R. Однако R предлагает значительно больше возможностей, чем мы показали здесь. Во второй части в этой серии из трех частей мы представим более продвинутые методы работы в R, начав с линейных и нелинейных регрессий.

Ресурсы

  • Примите участие в обсуждении материала на форуме.
  • Оригинал статьи (EN).
  • В следующей статье (EN) этой серии Брэд и Дэвид углубятся в функциональность R, используя возможности языка R для поиска и анализа аномальных данных.
  • Посетите главную страницу Проекта R для статистических вычислений (EN). Этот Web-сайт по R содержит обширную документацию, полные руководства и описания API. Для читателей, впервые встретившихся с R, особенно интересными будут следующие два документа:
  • Некоторые читатели колонки Дэвида Charming Python, будучи пользователями Python, высказывали особый интерес к привязкам Python к R. Фактически существует две такие привязки: RPy и более старая RSPython, которая тоже хороша (но Дэвид считает, что привязка RPy лучше). Каждая из них позволяет использовать весь диапазон функций R в коде на Python, используя объекты Python в качестве аргументов в функциях.
  • Дэвид написал в колонке Charming Python статью о Numerical Python (EN), который имеет похожий вид и многие из возможностей R (хотя R значительно обширнее). Вы можете посетить главную страницу Numerical Python..
  • Обратите внимание: на большинстве систем вы можете запустить браузер со сгенерированными HTML-страницами документации R, выполнив команду help.start() в командной строке R.
  • A Полная история S и S-Plus в Интернете.
  • Статья Камерона Лейрда Server clinic: R handy for crunching data (EN) (developerWorks, 2003) содержит хороший обзор R (равно как и S), а также список ресурсов.
  • Gnuplot может помочь вам с графиками и анализом данных. В статье Visualize your data with gnuplot (EN) (developerWorks, 2004), Нишан Састри представляет введение в Gnuplot, а также обсуждает проблему отсутствующих точек данных.
  • Стандарты для данных существуют, но мало кто о них знает. Например, если бы Брэд был разумнее, он мог бы с самого начала использовать формат данных Национального космического научного центра Common Data Format (CDF) Но это распространенная ошибка: если вы оказываетесь в подобной ситуации и вам нужно преобразовать исходные данные в определенный формат, подходящий для использования в средствах наподобие R, NumPy или Gnuplot, используйте такие средства, как Perl, Python или текстовые утилиты GNU. Дэвид написал руководство по последнему: Using the GNU text utilities (EN) (developerWorks, 2004).
  • Можете попрактиковаться в преобразовании и анализе данных о температуре, использованных Дэвидом и Брэдом в этой статье (скрипты на Python для конвертации исходного формата данных в более красивый табличный формат также здесь):(EN)
  • Много других ресурсов для Linux-разработчиков можно найти в разделе Linux сайта developerWorks.
  • Загрузите пробные версии ПО от IBM для Linux, включая WebSphere® Studio Application Developer, WebSphere Application Server, DB2® Universal Database, Tivoli® Access Manager и Tivoli Directory Server, воспользуйтесь руководствами и технической поддержкой раздела Speed-start your Linux app на developerWorks.(EN)
  • Присоединяйтесь к сообществу developerWorks, участвуя в блогах developerWorks.(EN)
  • Подберите книги по этим и другим техническим вопросам.(EN)

Комментарии

developerWorks: Войти

Обязательные поля отмечены звездочкой (*).


Нужен IBM ID?
Забыли Ваш IBM ID?


Забыли Ваш пароль?
Изменить пароль

Нажимая Отправить, Вы принимаете Условия использования developerWorks.

 


Профиль создается, когда вы первый раз заходите в developerWorks. Информация в вашем профиле (имя, страна / регион, название компании) отображается для всех пользователей и будет сопровождать любой опубликованный вами контент пока вы специально не укажите скрыть название вашей компании. Вы можете обновить ваш IBM аккаунт в любое время.

Вся введенная информация защищена.

Выберите имя, которое будет отображаться на экране



При первом входе в developerWorks для Вас будет создан профиль и Вам нужно будет выбрать Отображаемое имя. Оно будет выводиться рядом с контентом, опубликованным Вами в developerWorks.

Отображаемое имя должно иметь длину от 3 символов до 31 символа. Ваше Имя в системе должно быть уникальным. В качестве имени по соображениям приватности нельзя использовать контактный e-mail.

Обязательные поля отмечены звездочкой (*).

(Отображаемое имя должно иметь длину от 3 символов до 31 символа.)

Нажимая Отправить, Вы принимаете Условия использования developerWorks.

 


Вся введенная информация защищена.


static.content.url=http://www.ibm.com/developerworks/js/artrating/
SITE_ID=40
Zone=Linux, Open source
ArticleID=323396
ArticleTitle=Статистическое программирование на R: Часть 1. Купаемся в изобилии статистических возможностей
publish-date=07242008