Язык R. Генерация нормальной выборки из равномерной

[shukshin]

Есть функция, которая занимается генерированием нормальной выборки из равномерной.
(В учебных целях. Про 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) 
6. } 

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

PS lb <-> leftbound of uniform sample etc.
PPS про, вообще говоря, разное значение n в разных частях кода автор осведомлен
можно поставить $N = \operatorname{const}$ и звать

runif(N*n, lb, rb)

и так далее
Вопроса про эффективную реализацию это не снимает.

[provincialka]

shukshin
Не могли бы вы точнее объяснить, в чем проблема? Как я понимаю, вы реализуете численно ЦПТ? И что? медленно работает? Не то выдаёт?

-- 09.03.2016, 00:24 --

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

[_Ivana]

А разве в R нет библиотечной оптимизированной функции нормальнораспределенного рандома? В Матлабе и то есть.

[provincialka]

lb<-0;rb<-1
rn<-generate_norm(100,1,2)
hist(rn)

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

-- 09.03.2016, 00:34 --

_Ivana
ТС же написал:

shukshin в сообщении #1105159 писал(а):

В учебных целях. Про rnorm() все знают

[shukshin]

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$ элементов можно дождаться.

[Geen]

Может быть лучше другой алгоритм использовать?

[provincialka]

shukshin
На самом деле, не обязательно брать $N=100$, сумма равномерных вполне похожа на нормальную уже при $N=30$, работать будет в 3 раза быстрее.
И ещё. Зачем вы берёте произвольные lb,rb, а потом на каждом шаге используете их в вычислениях? Возьмите конкретные, например lb=-1,rb=1.

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

[provincialka]

Сделала пару модификаций функции (для 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

Принципиальной разницы нет. Но последняя функция не требует выделения дополнительной переменной (впрочем, она не глобальная, так что ничего страшного, если и создать)