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

в 16:16, , рубрики: Без рубрики

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

К мэйнстриму АйТи отнести себя не могу, но и «кнопки жму», если считать от первого школьного факультатива на Минск-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 раз вызываю по одному runif в каждом цикле for. Поверьте, это существенно для R.

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

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

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

Автор: vladob

Источник

Поделиться

* - обязательные к заполнению поля