GNU R: Вопросы и ответы

[Hogfather]

1. Трансформация данных

Возник вопрос: как преобразовать данные в R на манер сводных таблиц Excel. Вопрос был решен, но в процессе открыты для себя два полезных пакета, о которых хотелось бы рассказать.

Пакет reshape2 - позволяет агрегировать данные на манер Excel. На сайте http://www.r-statistics.com/tag/aggregate/ есть неплохой разбор его возможностей. Основная идея должна быть понятна вот из этой схемы (кликабельна).



Второй пакет -- sqldf. Позволяет писать прямо в коде SQL запросы к данным в синтаксисе SQLite.

Привожу простой код, который решает одну и ту же задачу с использованием двух этих пакетов. Имеются данные по весу кошачьих сердец, кошачьих тушек и пола (Пакет MASS данные cats), Попробуем найти для каждого пола число измерений и средний вес тушки.

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

Возможное решение.

> library(MASS)
> library(sqldf)
> library(reshape2)
> 
> head(cats)
  Sex Bwt Hwt
1   F 2.0 7.0
2   F 2.0 7.4
3   F 2.0 9.5
4   F 2.1 7.2
5   F 2.1 7.3
6   F 2.1 7.6
> 
> # Первый тест тест. Библиотека sqldf
> x<-Sys.time()
> for(i in 1:1000) z<-sqldf("select Sex,count(*) as cnt,avg(Bwt) as Avg_Bwt from cats group by Sex")
> y<-Sys.time()
> y-x
Time difference of 19.61596 secs
> z
  Sex cnt  Avg_Bwt
1   F  47 2.359574
2   M  97 2.900000
> # Второй тест. Библиотека reshape 2
> id<-1:nrow(cats)
> mcats<-melt(cbind(id,cats),id=c("id","Sex"))
> x<-Sys.time()
> for(i in 1:1000) {
+   df1<-dcast(subset(mcats,variable == "Bwt"),formula=Sex~variable,length)
+   df2<-dcast(subset(mcats,variable == "Bwt"),formula=Sex~variable,mean)
+   merge(df1,df2,intersect="Sex")
+ }
> y<-Sys.time()
> y-x
Time difference of 12.49125 secs
> z
  Sex cnt  Avg_Bwt
1   F  47 2.359574
2   M  97 2.900000

Как видно, второй вариант работает гораздо быстрее, но первый проще писать. Так что однозначного выбора нет.

Во втором примере также показано объединение двух таблиц по ключевому полю с использованием команды merge, поскольку команда cast не позволяет агрегировать сразу по двум функциям (может и позволяет, но я не умею).
Такие дела (с).


2. Графы в GNU R

Рассматривая интересные библиотеки в R можно упомянуть о возможности строить графы с помощью библиотеки igraph.

Простой пример кода. Пытаемся построить граф связей на портале аспирантов (фрагмент).

library(graph)
# Заполнем фрейм парами данных
left<-rep("Hogfather",16)
right<-c("Alextiger","caty-zharr","Dukar","Ink","IvanSpbRu","Martusya","osmos","saovu","Seta","Uzanka","Vica3","Димитриадис","Дмитрий В.","Домохозяйка","море","Степан Капуста")
left1<-rep("Степан Капуста",2)
right1<-c("Hogfather","Ink")
left<-c(left,left1)
right<-c(right,right1)
left1<-rep("море",10)
right1<-c("agasfer","Alextiger","bugo","fazotron","Hogfather","Ink","IvanSpbRu","Maksimus","Martusya","osmos")
left<-c(left,left1)
right<-c(right,right1)
myData<-data.frame(left,right)

# Формируем и строим граф
g<-graph.data.frame(myData,directed=F)
plot(g,vertex.label.cex=0.8,vertex.label.dist=1,vertex.size=10)

Результат



Нужно или нет, решать вам, но если надо рисовать серьезные графы, то рекомендую сперва глянуть в сторону graphviz. Очень кошерная и мощная штука.

[Hogfather]

Вышла книга. Шитиков В.К., Розенберг Г.С. Рандомизация и бутстреп: статистический анализ в биологии и экологии с использованием R. - Тольятти: «Кассандра», 2013. - 289 с.

Ссылка ведет на Интернет-версию на сайте авторов.

[Hogfather]

Сообщение от Diesel

Добрый день!
Опишу в кратце иеющиеся данные. В моём распоряжении есть пара тысяч текстовых файлов, в каждом из которых десяток тысяч строк вида: 1.2.2013 12:12 1234 1234 234 23423 (т.е. дата, время, четыре числа). В файлах содержится статистическая информация о работе некого сложного технического устройства. Один файл - это статистика работы одного устройства за период в несколько лет с интервалом в несколько минут.
Сейчас я провёл анализ данных десятка файлов. Определил несколько интересных закономерностей, для описания которых создал мат.модель. Однако в ручную обработать более 2 000 файлов - процесс долгий, т.е. нужна автоматизация. Плюс есть множество идей по выявлению других закономерностей Анализ нужно проводить как по каждому отдельному файлу, т.к. и по группе файлов.

Не рассматривая собственно вопрос анализа, продемонстрирую как можно закачать данные в пакет R. Для тестового примера был создан каталог Data_R, в которые положены 2 одинаковых файла test1.txt и test2.txt, вида:

1.2.2013 12:12 1234 1234 234 23423
1.2.2013 12:12 1234 1234 234 23423
1.2.2013 12:12 1234 1234 234 23423
1.2.2013 12:12 1234 1234 234 23423
1.2.2013 12:12 1234 1234 234 23423
1.2.2013 12:12 1234 1234 234 23423

Простенький скрипт, который создает список всех txt файлов в данном каталоге, а потом импортирует их во внутреннюю таблицу R, добавляя колонку "Источник" прилагаю.

FileList<-dir("Data_R",pattern = "*.txt", full.names = TRUE, ignore.case = TRUE, include.dirs =TRUE)
MyResult<-data.frame()
for(FileName in FileList)
{
  TempTable<-read.table(FileName,sep=" ")
  TempTable$FileName<-FileName
  MyResult<-rbind(MyResult,TempTable)
}

Для просмотре результатов работы вводит команду MyResult, которая, собственно, отобразит все значения этой таблицы

> MyResult
V1 V2 V3 V4 V5 V6 FileName
1 1.2.2013 12:12 1234 1234 234 23423 Data_R/test1.txt
2 1.2.2013 12:12 1234 1234 234 23423 Data_R/test1.txt
3 1.2.2013 12:12 1234 1234 234 23423 Data_R/test1.txt
4 1.2.2013 12:12 1234 1234 234 23423 Data_R/test1.txt
5 1.2.2013 12:12 1234 1234 234 23423 Data_R/test1.txt
6 1.2.2013 12:12 1234 1234 234 23423 Data_R/test1.txt
7 1.2.2013 12:12 1234 1234 234 23423 Data_R/test2.txt
8 1.2.2013 12:12 1234 1234 234 23423 Data_R/test2.txt
9 1.2.2013 12:12 1234 1234 234 23423 Data_R/test2.txt
10 1.2.2013 12:12 1234 1234 234 23423 Data_R/test2.txt
11 1.2.2013 12:12 1234 1234 234 23423 Data_R/test2.txt
12 1.2.2013 12:12 1234 1234 234 23423 Data_R/test2.txt

Как видим, все работает. Дальше можно уже делать срезы какие душе угодно.

Добавлено через 36 минут
С другой стороны, данных много, поэтому возможно удобнее будет загнать в базу данных SQLite (см. выше о работе в R) и использовать срезы оттуда.

Модифицируем чуть-чуть вышеприведенный код

require(RSQLite)
drv <- dbDriver("SQLite")
fname<-"mytest.sqlite"
conn <- dbConnect(drv, dbname = fname)
if(dbExistsTable(conn, "RESULTS")) {
dbRemoveTable(conn, "RESULTS")
}
FileList<-dir("Data_R",pattern = "*.txt", full.names = TRUE, ignore.case = TRUE, include.dirs =TRUE)
for(FileName in FileList)
{
TempTable<-read.table(FileName,sep=" ")
TempTable$FileName<-FileName
if(dbExistsTable(conn, "RESULTS")) {
dbWriteTable(conn, "RESULTS", TempTable, append = T)
}
else
{
dbWriteTable(conn, "RESULTS", TempTable)
}
}
dbDisconnect(conn)
dbUnloadDriver(drv)

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



Добавлено через 49 минут

Сообщение от Kayra

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

Kayra, Пример кода.

set.seed(666)
age <- 1:10
y.low <- rnorm(length(age), 150, 25) + 10*age
y.high <- rnorm(length(age), 250, 25) + 10*age
plot(age,y.high,type = 'n', ylim = c(100, 400), ylab = 'Y Range', xlab = 'Age (years)')
lines(age, y.low )
lines(age, y.high)
polygon(c(age, rev(age)), c(y.high, rev(y.low)), col = "lightblue", border = NA)

# Второй график, чтобы поржать

y.low <- rnorm(length(age), 150, 25) + 10*age
y.high <- rnorm(length(age), 250, 25) + 10*age
lines(age, y.low, col = "green" )
lines(age, y.high, col = "green")
polygon(c(age, rev(age)), c(y.high, rev(y.low)), col = "green", border = NA)

[Uzanka]

Hogfather,
а в R можно оценивать GARCH модели с разными распределениями ошибок?

а оценивать Stochastic volatility models с разными распределениями?

Или самим код надо писать?

ЗЫ. А что-то типа фильтра Калмана там есть?

[Hogfather]

Сообщение от Uzanka

а в R можно оценивать GARCH модели с разными распределениями ошибок?

Для примера потребуются две библиотеки fGarch и evir
Пример из книжки, упомянутой ниже

library(fGarch)
data(bmw,package="evir")

bmw.garch_norm = garchFit(~arma(1,0)+garch(1,1),data=bmw,cond.dist="norm")
options(digits=3)
summary(bmw.garch_norm)
options(digits=10)


x = bmw.garch_norm@residuals / bmw.garch_norm@sigma.t
n=length(bmw)
grid = (1:n)/(n+1)
fitdistr(x,"t")

par(mfrow=c(1,2))
qqnorm(x,datax=T,ylab= "Standardized residual quantiles", 
       main="(a) normal plot",
       xlab="normal quantiles")
qqline(x,datax=T)
qqplot(sort(x), qt(grid,df=4),
       main="(b) t plot, df=4",xlab= "Standardized residual quantiles",
       ylab="t-quantiles")
abline(   lm(   qt(c(.25,.75),df=4)~quantile(x,c(.25,.75))   )   )

bmw.garch_t = garchFit(~arma(1,1)+garch(1,1),cond.dist="std",data=bmw)
options(digits=4)
summary(bmw.garch_t)

Справочник по финансовым пакетам: http://cran.r-project.org/web/views/Finance.html

Рекомендую почитать книжку Statistics and Data Analysis for Financial Engineering. Я купил и мои волосы теперь чистые и шелковистые.

Сообщение от Uzanka

а оценивать Stochastic volatility models с разными распределениями?

http://cran.r-project.org/web/packag...vol/index.html

Сообщение от Uzanka

ЗЫ. А что-то типа фильтра Калмана там есть?

Пакет
http://cran.r-project.org/web/packages/FKF/index.html
Статья по теме
http://www.jstatsoft.org/v39/i02/paper

[Uzanka]

Hogfather,
просто огромное спасибо!!!!!!! пошла изучать

[Hogfather]

Приобрел неплохую книжку по Data Mining в GNU R, рекомендую.




Пакет rattle предназначен для поиска закономерностях в данных (Data Mining) с помощью регрессионных деревьев, кластерного анализа и метода опорных векторов. В книжке разбирается порядок действий от загрузки данных до интерпретации результатов.

library(rattle)
rattle()

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

UPD: Появилось сообщение на форуме поддержки. На настоящий момент решение это проблемы выглядит вот так [1] . Вложение со старой функцией удалено, чтобы не смущало.

The pmmlCanExport is an issue with the new release of the pmml package - it is being addressed at present - a quick fix:

> pmmlCanExport <- function(...) TRUE
> rattle()

Regards,

Graham Williams
http://togaware.com

Т.е. есть еще некоторые проблемки с продуктом,а так, в целом, здорово!

[Hogfather]

UPD: Под МакОс rattle категорически отказался работать. Проблемы с GTK+

[Hogfather]

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

Сообщение от Paul Kellerman

Имеет тест на 45 вопросов по 4 варианта ответа в каждом. Только один вариант верный.
Тест разбит на 3 блока по 15 вопросов. За правильные ответы на вопросы из 1-го блока
начисляется 1 балл, из 2-го блока - 2 балла, из 3-го блока - 3 балла, и соответственно,
максимальный балл, который можно набрать = 90. Определить вероятность сдачи теста,
тыкая ответы наугад, при условии, что для успешной сдачи нужна набрать минимум 61.

А давайте её еще немного усложним и ко всему прочему построим графики функции вероятности и функции распределения?

[Paul Kellerman]

Сообщение от Paul Kellerman

Задан тест на 45 вопросов по 4 варианта ответа в каждом. Только один вариант верный.
Тест разбит на 3 блока по 15 вопросов. За правильные ответы на вопросы из 1-го блока
начисляется 1 балл, из 2-го блока - 2 балла, из 3-го блока - 3 балла, и соответственно,
максимальный балл, который можно набрать = 90. Определить вероятность сдачи теста,
тыкая ответы наугад, при условии, что для успешной сдачи нужна набрать минимум 61.

Сообщение от Hogfather

А давайте её еще немного усложним и ко всему прочему построим
графики функции вероятности и функции распределения?

Пусть T - дискретная случайная величина, равная количеству набранных баллов (0...90).
Функция вероятности случайной величины T: f(t) = P(T=t)
Функция распределения случайной величины T: F(t) = P(T<=t)

Решение в математической среде Waterloo Maple 15.0.