В Хаброжители я прибился относительно недавно.
К мэйнстриму АйТи отнести себя не могу, но и «кнопки жму», если считать от первого школьного факультатива на Минск-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
