Практическое приложение к «I Love R». Неплановое

В Хаброжители я прибился относительно недавно.

К мэйнстриму АйТи отнести себя не могу, но и «кнопки жму», если считать от первого школьного факультатива на Минск-22, уже, почитай,… 40 лет… (О, боже! столько же не живут!)

Ах да, я отвлекся…
Написал я дрожащими (от волнения) пальцами свой первый пост про R и был весьма польщен полученными отзывами. Что еще более подкрепило желание выполнить обещание и показать кое что из практического применения R. И в частности из области биоинформатики, где R наиболее популярен, и где мы вместе с ним и трудимся.

В то же время, уже не впервые вижу, что R применяют в качестве языка для «макетирования», с целью потом переписать на чем-нибудь «настоящем» — С++, ну или, на худой конец, на Python’е.

Вот конкретно этот пост был своего рода спровоцирован статьей про индексный метод вычисления вероятностей. Оставляю сейчас свои придирки к изложению в статье касаемо теории вероятностей (что тоже непросто, имея за плечами 20+ лет преподавания, в т.ч. тервера).
Под катом — пример приведения R-кода из этой статьи к «рабочему» (по моему субъективному мнению) состоянию.

Вкратце:
Уважаемый kokorins написал полезную по моему мнению статью про индексный способ построения генератора дискретной случайной величины.
В статье был использован вот такой код, чтобы проиллюстрировать «наивную» версию алгоритма:

#Naive approach len<-10 ps<-seq(len,2, by=-1) ps<- 1/ps^2 ps<-ps/sum(ps) ss<-cumsum(ps) gen_naiv <- function() {   alpha<-runif(1)   return (min(which(alpha<ss))) }  #sample val<-c() len<-10000 for(i in 1:len) {   val<-append(val, gen_naiv()) } vals<-factor(val) plot(vals) 

Меня заинтересовал алгоритм, но ключевые слова, которые я «выудил» из цитируемого поста оказались малоспецифичны, а фамилия автора — Chen — вообще чуть ли не самая распространенная в Китае, а посему одна их самых многочисленных в интернете…

Авторский «стиль» программирования в R далек от совершенства (мягко говоря), но статья была не про R. Я сначала пытался не обращать внимание. Но в попытках разобраться с алгоритмом по авторскому R-тексту я потерял терпение и решил использовать этот текст в качестве показательного примера.

Итак, как уже сказано — код рабочий. Но очень плохой. И я попробую показать чем, почему, и как сделать лучше.

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

len<-10 ps<-seq(len,2, by=-1) ps<- 1/ps^2 ps<-ps/sum(ps) ss<-cumsum(ps)  gen_naiv <- function() {   alpha<-runif(1)   return (min(which(alpha<ss))) }  #sample get_naive<-function(len=10000){ 	val<-c() 	st<-system.time(  	for(i in 1:len) {   		val<-append(val, gen_naiv()) 	}) 	print(st) 	return(val) } 

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

Выполнив в среде R указанный код, я теперь могу получить последовательности разной длины, заодно осведомляясь о затраченном на выполнение времени.

> get_naive(10)    user  system elapsed    0.006   0.000   0.007   [1] 8 8 8 8 9 5 8 9 9 7 > get_naive(100)    user  system elapsed    0.008   0.000   0.008    [1] 8 1 7 8 8 4 9 8 7 8 9 9 9 9 5 7 9 9 8 3 7 1 1 9 7 7 9 9 5 8 6 6 4 9 9 6 9 8 9 8 6 8 6 9 2 8 9 9 9 9 7 8 5 9 6 6 5 9 9 8 9 9 7 8 6 8 9 1 8 9  [71] 9 3 5 9 8 9 9 9 9 6 9 9 5 9 8 7 8 9 9 9 8 9 7 9 8 7 8 5 5 8 >   

Полюбовавшись на тысячные доли секунды, затраченные на выполнение и на собственно псевдослучайную последовательность чисел от 1 до 9, «спрячу» вывод функции в переменнную val, но продолжу эксперименты по наращиванию длины сэмпла:

> val<-get_naive(1000)    user  system elapsed    0.021   0.000   0.021  > val<-get_naive(10000)    user  system elapsed    0.342   0.011   0.354  > val<-get_naive(100000)    user  system elapsed    24.86    8.55   33.48  

Ай-ай-ай!
Похоже, непорядок! Как-то уж очень круто растет время выполнения.
Может это действительно проблемы наивного подхода? А алгоритм имени китайского товарища работает в разЫ быстрее?

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

> get_chen(1000)    user  system elapsed    0.031   0.001   0.032  > get_chen(10000)    user  system elapsed    0.508   0.063   0.567  > get_chen(100000)    user  system elapsed   46.303   5.946  53.088  

Оп-па! Вот так неожиданность… Подкачал, товарищ Чен. Что-то «улучшенный» алгоритм отказывается работать лучше «наивного». Или это цитируемый мной автор не смог донести весь «блеск»? Я специально не стал показывать очередную порцию плохого по моему мнению кода. Лучше не привыкать.

Мой первоначальный план был шаг за шагом разобрать ошибки коллеги и показать как надо. Но я решил просто переписать с комментариями код («наивную» версию), как бы это сделал я, опираясь на свой опыт работы с R и оставить kokorins‘ a самого бороться со своими недостатками или недостатками своего или китайского кода. Извините меня — просто не смог.

Вот код, который у меня получился:

M<-10             # Как у автора ps<- 1/(M:2)^2    # вектор вероятностей появления отдельных значений; пока в каких-то относительных единицах ps<-ps/sum(ps)    # Нормируем сумму вероятностей на единицу ss<-c(0,cumsum(ps))    # Строим вектор границ инервалов  lbs=1:(M-1)  # Это "метки на кубике" (вектор) -- те самые дискретные значения, которые пойдут на выход   get_vladob <- function(n, brks=ss, labs=lbs) { 	st<-system.time(   		val<-as.numeric(cut(runif(n), breaks=brks, labels=labs))   	)   	print(st)   	return(val) } 

И результаты «прогона»

> val<-get_vladob(1000)    user  system elapsed    0.010   0.000   0.012  > val<-get_vladob(10000)    user  system elapsed    0.004   0.000   0.004  > val<-get_vladob(100000)    user  system elapsed    0.036   0.001   0.038  > val<-get_vladob(1000000)    user  system elapsed    0.581   0.029   0.606  > val<-get_vladob(10000000)    user  system elapsed    4.015   0.311   4.318  > val<-get_vladob(100000000)    user  system elapsed   38.405   3.092  41.578  

Как видите — совсем неплохо для «медленного» языка.

Я сейчас с интересом посматриваю в сторону других языков. С интересом увидел бы результаты подобного теста, реализованного в других средах. В частности, мне жутко любопытно получит ли коллега kokorins результаты на С++ и какие.

Что поменялось в коде?
В первую очередь, я использовал «векторные» свойства R. Например, если мне нужны 1000 значений функции runif, я использую runif(1000) и векторные выражения для обработки, а не 1000 раз вызываю по одному funif в каждом цикле for. Поверьте, это существенно для R.

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

В третьих, я упоминал в своем первом посте, что R написан «статистиками для статистиков». Функция cut назначает метки (labels) значениям в векторе в зависимомости от того, в какой интервал эти значения попадают. Вектор границ интервалов передается в функцию cut, как параметр breaks. Это достаточно распространенная операция при обработке данных (группировка, например), чтобы быть реализованной (и оптимизированной!) в R в виде отдельной функции.

В общем — R.
Прошу любить и жаловать