2014 dxdy logo

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

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




Начать новую тему Ответить на тему
 
 Язык R. Генерация нормальной выборки из равномерной
Сообщение08.03.2016, 23:22 


20/10/12
235
Есть функция, которая занимается генерированием нормальной выборки из равномерной.
(В учебных целях. Про 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 
Заслуженный участник
Аватара пользователя


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

-- 09.03.2016, 00:24 --

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

 Профиль  
                  
 
 Re: Знатокам языка R
Сообщение09.03.2016, 00:31 


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

 Профиль  
                  
 
 Re: Знатокам языка R
Сообщение09.03.2016, 00:33 
Заслуженный участник
Аватара пользователя


18/01/13
12044
Казань
Код:
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 


20/10/12
235
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 
Заслуженный участник
Аватара пользователя


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

 Профиль  
                  
 
 Re: Язык R. Генерация нормальной выборки из равномерной
Сообщение09.03.2016, 14:21 
Заслуженный участник
Аватара пользователя


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

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

 Профиль  
                  
 
 Re: Язык R. Генерация нормальной выборки из равномерной
Сообщение09.03.2016, 15:27 
Заслуженный участник
Аватара пользователя


18/01/13
12044
Казань
Сделала пару модификаций функции (для 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 ] 

Модераторы: Karan, Toucan, PAV, maxal, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group