2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




 
 Язык R. Генерация нормальной выборки из равномерной
Сообщение08.03.2016, 23:22 
Есть функция, которая занимается генерированием нормальной выборки из равномерной.
(В учебных целях. Про rnorm() все знают)
  1. generate_norm <- function(n, mu, sigma){ 
  2.   uni_matrix <- matrix(runif(n^2, lb, rb), n) 
  3.   out <- apply(uni_matrix, 1, sum) 
  4.   out <- sigma*((out - n * ((lb + rb)/2.0)) / sqrt(n * ((rb - lb)^2/12.0))) + mu 
  5.   return(out) 

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

PS lb <-> leftbound of uniform sample etc.
PPS про, вообще говоря, разное значение n в разных частях кода автор осведомлен
можно поставить $N = \operatorname{const}$ и звать
Код:
runif(N*n, lb, rb)
и так далее
Вопроса про эффективную реализацию это не снимает.

 
 
 
 Re: Знатокам языка R
Сообщение09.03.2016, 00:22 
Аватара пользователя
shukshin
Не могли бы вы точнее объяснить, в чем проблема? Как я понимаю, вы реализуете численно ЦПТ? И что? медленно работает? Не то выдаёт?

-- 09.03.2016, 00:24 --

Я запустила вашу функцию. Только пришлось заранее задать параметры lb, rb. Какой-то ответ даёт.
Кстати, если "лагает" скрипт, то, может, дело в нём, а не в функции?

 
 
 
 Re: Знатокам языка R
Сообщение09.03.2016, 00:31 
А разве в R нет библиотечной оптимизированной функции нормальнораспределенного рандома? В Матлабе и то есть.

 
 
 
 Re: Знатокам языка R
Сообщение09.03.2016, 00:33 
Аватара пользователя
Код:
lb<-0;rb<-1
rn<-generate_norm(100,1,2)
hist(rn)

Гистограмма выглядит вполне "нормально". :D

-- 09.03.2016, 00:34 --

_Ivana
ТС же написал:
shukshin в сообщении #1105159 писал(а):
В учебных целях. Про rnorm() все знают

 
 
 
 Re: Знатокам языка R
Сообщение09.03.2016, 10:39 
provincialka
да, медленно работает для выборок > 10000. (я просто не знаю насколько быстр R, может это и нормально)
учитывая, что я пишу на R третий день, я предположил, что можно написать лучшую реализацию

а я вот так сделал и впринципе стало лучше (так даже логичнее стало) :
Код:
generate_norm <- function(n, mu, sigma){
  N = 100
  uni_matrix <- matrix(runif(N*n, lb, rb), n)
  out <- apply(uni_matrix, 1, sum)
  out <- sigma*((out - N * ((lb + rb)/2.0)) / sqrt(N * ((rb - lb)^2/12.0))) + mu
  return(out)
}


уже и $10^5$ элементов можно дождаться.

 
 
 
 Re: Знатокам языка R
Сообщение09.03.2016, 11:23 
Аватара пользователя
Может быть лучше другой алгоритм использовать?

 
 
 
 Re: Язык R. Генерация нормальной выборки из равномерной
Сообщение09.03.2016, 14:21 
Аватара пользователя
shukshin
На самом деле, не обязательно брать $N=100$, сумма равномерных вполне похожа на нормальную уже при $N=30$, работать будет в 3 раза быстрее.
И ещё. Зачем вы берёте произвольные lb,rb, а потом на каждом шаге используете их в вычислениях? Возьмите конкретные, например lb=-1,rb=1.

Но всё-таки уточните задачу. Надо построить "нормально распределённую величину" именно на основе ЦПТ? Или просто через равномерное распределение?

 
 
 
 Re: Язык R. Генерация нормальной выборки из равномерной
Сообщение09.03.2016, 15:27 
Аватара пользователя
Сделала пару модификаций функции (для N=30).
код: [ скачать ] [ спрятать ]
Используется синтаксис Python
generate_norm1 <- function(n, mu, sigma){
  uni_matrix <- matrix(runif(30*n, -1, 1), n)
  out <- apply(uni_matrix, 1, sum)
  out <- out*sigma/sqrt(10) + mu
  return(out)
}

generate_norm2 <- function(n, mu, sigma){
  uni_matrix <- runif(30*n, -1, 1)
  dim(uni_matrix) <- c(n,30)
  out <- apply(out, 1, sum)
  out <- out*sigma/sqrt(10) + mu
  return(out)
}

generate_norm3 <- function(n, mu, sigma){
  out <- runif(30*n, -1, 1)
  dim(out) <- c(n,30)
  out <- apply(out, 1, sum)*sigma/sqrt(10) + mu
  return(out)
}

print(system.time(rn<-generate_norm1(n,1,2)))
print(system.time(rn<-generate_norm2(n,1,2)))
print(system.time(rn<-generate_norm3(n,1,2)))


Вот результаты по времени просчёта

При n=1000000
    пользователь система прошло
    11.78 0.31 12.31
    пользователь система прошло
    10.85 0.25 11.33
    пользователь система прошло
    10.92 0.23 11.41
Принципиальной разницы нет. Но последняя функция не требует выделения дополнительной переменной (впрочем, она не глобальная, так что ничего страшного, если и создать)

 
 
 [ Сообщений: 8 ] 


Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group