Портал аспирантов - 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

Цитата:

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

Итак, вернемся к нашим баранам. Собственно, решение задачи основано на знании формулы Бернулли и понимания биноминального распределения.

Функция вероятности для биноминального распределения в R имеет вид dbinom(x, size, prob, log = FALSE), остальные связанные функции:
pbinom(q, size, prob, lower.tail = TRUE, log.p = FALSE) - функция распределения
qbinom(p, size, prob, lower.tail = TRUE, log.p = FALSE) - квантили распределения
rbinom(n, size, prob) - генерация случайного числа, распределенного по данному закону
где, соответственно: size - число попыток, prob - вероятность удачного исхода

Поэтому задачу:

Цитата:

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

Вот к примеру есть тупо тест тупо на 45 вопросов тупо по 4 ответа
в каждом, и только один ответ верный. Для положительной оценки
нужно минимум на 31 вопрос ответить верно. Какова вероятность
сдать его, ничего не зная, тупо тыкая наугад ответы?

можно решить двумя способами, через сложение устраивающих нас исходов испытания или как разницу между 1 и вероятностью неустраивающих нас результатов, которая равна куммулятивной сумме вероятностей исходов от 0 до 30, что соответствует значению функции распределения pbinom(30, 45, 0.25)

Код:

> sum(dbinom(31:45, 45, 0.25))

[1] 7.527228e-10

> 1-pbinom(30, 45, 0.25)

[1] 7.527228e-10

В MS Excel это считается так:
=1-БИНОМ.РАСП(30;45;0,25;ИСТИНА)
и имеем результат
7,52723E-10

Вооружившись этими немудреными знаниями приступаем к задаче. Здесь у нас существуют возможные варианты исходов от 0 до 90, но поскольку веса у каждого блока в тесте разные, проще прогнать все варианты в цикле и учесть число возможных перестановок, а также вероятность совместного события. Поскольку в данном случае события независимые, то вероятности перемножаем. По сути, мы решаем ту же задачу, только слегка усложненную.

Код:

> # Инициируем переменные

> P<-0   # Переменная для подсчета куммулятивной суммы, как увидим дальше, она лишняя. Просто для явности подсчетов для решения задачи

> Z<-data.frame(X=0:90,y=rep(0,91)) # Здесь мы будем накапливать вероятности для исходо теста 0:90

> for(i in 0:15) for (j in 0:15) for(k in 0:15) { # Перебираем возможные варианты для каждой части теста

+   Z$y[i+2*j+3*k+1]<-Z$y[i+2*j+3*k+1]+dbinom(i,15,0.25)*dbinom(j,15,0.25)*dbinom(k,15,0.25)  # считаем вероятность куммулятивно для точки

+   if(i+2*j+3*k>=61) P<-P+dbinom(i,15,0.25)*dbinom(j,15,0.25)*dbinom(k,15,0.25) # считаем для точки

+ }

> Z$cdf<-cumsum(Z$y)

> cat("Ответ на задачу составляет: ",P)

Ответ на задачу составляет:  1.73944e-08

> oldpar<-par(mfrow=c(2,1))

> plot(Z$y~Z$X,main="Функция вероятности f(x)",xlab="x",ylab="f(x)",pch=19,cex=0.8,col="blue")

> points(Z$y[62:91]~Z$X[62:91],cex=0.9,col="red",pch=19)

> text(70,0.02,expression(sum(f[i],i==61,90)==1.73944 %*% 10^{-8}))

> plot(Z$cdf~Z$X,main="Рапределения F(x)",xlab="x",ylab="f(x)",pch=19,cex=0.8,col="blue")

> abline(h=Z$cdf[61],v=60,col="darkgrey",lty=2)

> points(Z$cdf[61]~Z$X[61],cex=0.9,col="red",pch=19)

> par(oldpar)

> # Собственно, как и было сказано выше, вероятности также можно посчитать как:

> 1-Z$cdf[61]  # 1-F(60); 61, а не 60, потому как нумерация в массиве начинается с 1

[1] 1.73944e-08

> sum(Z$y[62:91]) # sum(f(t),t = 61...90)

[1] 1.73944e-08

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

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

о так вот вы куда ушли с решением !
пик и есть нужный нам процент?

Цитата:

Сообщение от _Tatyana_ (Сообщение 404329)

пик и есть нужный нам процент?

Раз есть вопросы, значит график неудачный. Перерисовал.

Итак, пик на верхнем графике соответствует 22 баллам со значением вероятности 6.315466e-02 или 6%. Т.е. случайно в 6% случаев можно ответить на 22 балла. Вообще, верхний график показывает вероятности набрать случайным образом баллы от 0 до 90. Поэтому, если нам надо посчитать вероятность не менее 61, мы складываем вероятности набрать 61,62, ... , 90 баллов (точки выделены красным) и получаем нужный нам ответ 1.73944e-08.

Нижний график показывает куммулятивную сумму вероятностей, т.е. в точке 60 (красная точка) он равен сумме вероятностей набрать баллы от 0 до 60, иными словами, значение в каждой точке показывает вероятность набрать не более заданного количества баллов. Для того, чтобы найти вероятность набрать более чем заданное количество баллов (в нашем случае 60, потому как 61>60 и по условиям задачи оно входит в вопрос), учитывая, что сумма всех вероятностей равна 1, имеем 1-F(60)=1.73944e-08

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

Цитата:

Сообщение от _Tatyana_ (Сообщение 404569)

а вопрос задачи был о 31 балле. распределение говорит что это три процента
правильно?

_Tatyana_, вроде бы там написано

Цитата:

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

1.73944e-08