Рубрика: R

Как запускать R скрипты из SAS?

Всем привет.

r_logo

В сегодняшнем уроке Вы узнаете как запускать R скрипты непосредственно из SAS.

Читать далее Как запускать R скрипты из SAS?

Чтение SAS data sets из R

Допустим Вы подготовили данные для анализа в SAS и хотите воспользоваться богатыми возможностями предоставляемыми многочисленными R-библиотеками. Причин на это может быть несколько.
  • Например у вас нет лицензии на нужный продукт (IML, Miner, и т.д.)
  • В SAS ещё нет соответствующего метода/алгоритма/функционала (удобный синтаксис для операций над матрицами, какие-то новые методы изобретенные недавно, графические возможности (ggplot2), парсинг различных форматов данных, подготовка отчетов и презентаций (knitr) и т.д,)
  • Вы хотите иметь доступ к source code для того что бы разобраться детально как работает алгоритм и, возможно, изменить его под свои нужды
Последний пункт, на мой взгляд, особенно важен. Ну и конечно ни один продукт/язык программирования не решает все задачи одинаково хорошо, поэтому иногда имеет смысл использовать сильные стороны имеющихся в распоряжении интсрументов. А R — интсрумент мало того что хороший, но еще и бесплатный.
Со своей стороны SAS предоставляет интерфейс для вызова R функций через IML:
Calling Functions in the R Language. Но его у Вас может не быть, поэтому можно экспортировать данные в CSV формат и потом считать их в R. Более быстрый и удобный способ — воспользоваться библиотекой haven. Недавно вышла новая версия этого пакета, о чем автор написал на сайте r-bloggers.
Пример чтения:
# Установка package:
install.packages("haven")
# Загрузка библиотеки в R
library(haven)
Описание функций пакета можно найти на сайте CRAN в pdf файле.
Пример входных данных можно скачать из GitHub репозитория: hadley.sas7bdat. Сохранив файл локально, считываем его командой:
data = read_sas("D:/_tmp/hadley.sas7bdat")
На выходе получаем data.frame структуру данных с именованными переменными:
​

Простой пример визуализации данных в R

Что-то давно никто не писал тут, походу все в отпусках. Раз так, то позволю себе написать пост про применение R для построения графиков. Данные будут взяты из предыдущего поста. Для этого можно пройти на страницу Wiki, скопировать и сохранить табличку в csv file: link.

Считываются данные одной строкой:


data <- read.table(«./AnscombeQuartet.csv», header = TRUE, sep =»,»)

После этого все столбцы доступны по соответствующему header name:

data$x1

Далее можно проверить статистики, соответствуют ли они заявленным значениям.
Сделать это можно двумя способами:

  1. Применить к каждому столбцу функцию: mean(data$x1), var(data$x1), ….
  2. Либо воспользоваться функцией sapply(). Это некий аналог функции map в функциональных языках программирования.
Мы воспользуемся вторым способом 
(sapply(data, mean); sapply(data, var)):

Все статистики совпадают с заявленными. Теперь построим графики:
Если нужно построить всего один график, то делается это очень просто;

plot(data$x1, data$y1, xlab = «x», ylab = «y») #xlab, ylab add label names

Но нам нужно пробежаться по всем столбцами и уместить все 4 графика в одну фигуру. Это можно сделать в цикле:

par(mfrow=c(2,2), mar=c(2.0, 1.5, 1.5, 1.5), oma = c(2,2,2,2), pch = 16)
for(i in 1:4){
  x = eval(parse(text = paste( «data$x», i, sep=»» )))
  y = eval(parse(text = paste( «data$y», i, sep=»» )))
  plot( x, y )
  abline( lm(y ~ x) )
  grid( 10,10 )
}

Получаем график (на всякий случай в высоком разрешении: URL):

Пояснения:

  • abline( lm(y ~ x) ) строит линейную модель. Сама функция lm() возвращает объект определенного класса, и для того чтобы построить выходные данные не нужно их специальным образом извлекать, потому как функция abline( ) сама распознает контекст. 
  • Функция grid() добавляет вертикальные и горизонтальные линии координатной сетки
  • Функцией par( ) задаются расстояния от границ каждого графика до границ всей фигуры, также расстояния между своими фигурами, pch — вид маркера. Но можно было бы использовать dafault значения.
  • mfrow() задает геометрию расположения графиков
Вообще в R очень и очень богатые графические возможности (помимо статистики и более 5000 packages). Даже поддерживаются несколько парадигм, например есть ggplot2. 

Из SAS R code можно запустить очень просто, подробное описание есть ТУТ.


Весь R код использованный при написании поста можно найти ТУТ.

Об интервалах значений функции плотности вероятности

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

 

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

[ f( x | mu, sigma^2 ) = frac{1}{sigma sqrt{2 pi}} exp(-frac{1}{2} (frac{x-mu}{sigma})^2 ) ]

Интуитивно я все время думал, что раз уж речь идет о вероятности, то значение ну никак не может быть больше 1. Изначально, первый пост я хотел написать про Байесовые классификаторы, но он еще не окончен. Так вот, вычисляя условные вероятности, во время имплементации классификатора, я столкнулся с этим явлением, и стал искать ошибку, ну а потом провел простой эксперимент.

Чтобы посчитать плотность вероятность в точке, запишем функцию:

probability <- function(x, mu, sigma) {
    p <- (1/(sigma * sqrt(2 * pi))) * exp(-0.5 * ((x - mu)/sigma)^2)
    return(p)
}
probability(5.1, 5.06, 0.35)
## [1] 1.132

То же самое можно сделать проще, используя встроенную функцию:

dnorm(5.1, 5.06, 0.35)
## [1] 1.132

Оказывается это совсем не ошибка, а хорошо известная вещь. Картина проясняется, если вспомнить, что речь идет именно о плотности веротности, которая действительно может принимать значения больше единицы. На просторах интернета я нашел следующий поясняющий пример. Рассмотрим непрерывное равноемерное распределение на интервале ( [0, frac{1}{2}] ) с плотностью вероятности в каждой точке интервала равным ( 2 ) и ( 0 ) во всех остальных точках. Тогда интеграл от ( 0 ) до ( 2 ) равен ( 1 ), как того и требует определение.