Портал аспирантов - GNU R: Вопросы и ответы

Портал аспирантов (http://www.aspirantura.spb.ru/forum/index.php)

- Software (программное обеспечение) (http://www.aspirantura.spb.ru/forum/forumdisplay.php?f=107)

- - GNU R: Вопросы и ответы (http://www.aspirantura.spb.ru/forum/showthread.php?t=10501)

Трансформация данных в R и графы

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

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

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

http://www.r-statistics.com/wp-conte...st-300x214.png

Второй пакет -- 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)

Результат

http://www.aspirantura.spb.ru/forum/...pictureid=1100

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

Цитата:

Сообщение от Diesel (Сообщение 356476)

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

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

http://aspirantura.spb.ru/forum/pict...pictureid=1162

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

Цитата:

Сообщение от Kayra (Сообщение 356490)

В какой программе можно построить такие закрашенные полосы?
Пара кривых строится по точкам 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)

http://aspirantura.spb.ru/forum/pict...pictureid=1163

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

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

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

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

Цитата:

Сообщение от Uzanka (Сообщение 396699)

а в 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 (Сообщение 396699)

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

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

Цитата:

Сообщение от Uzanka (Сообщение 396699)

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

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

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

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

http://ecx.images-amazon.com/images/...SH20_OU01_.jpg

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

Код:

library(rattle)

rattle()

http://aspirantura.spb.ru/forum/pict...pictureid=1364

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

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

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

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

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

Цитата:

Сообщение от Paul Kellerman (Сообщение 403562)

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

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

Цитата:

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

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

Цитата:

Сообщение от Hogfather (Сообщение 403606)

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

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

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

http://4put.ru/pictures/max/772/2372608.jpg