Быстрая генерация 1000 средних точек выборки из усеченного гамма-распределения с 1000 различными значениями форм и масштабов в R
После поиска по форуму я не нашёл похожих вопросов. Если вы найдете один, пожалуйста, дайте мне знать. Я был бы очень признателен.
Мне нужно сгенерировать 1000 средних точек выборки из усеченного гамма-распределения с 1000 различными значениями форм и масштабов в R.
Мой следующий код работает, но очень медленно. Как улучшить производительность?
library(distr)
library(distrEx)
library(truncdist)
set.seed(RANDOM.SEED)
shape.list <- runif(1000, max = 10, min = 0.01)
scale.list <- runif(1000, max = 100000, min = 100000)
mean.list <- list()
std.dev.list <- list()
for (i in seq(1000)) # very slow
{
sample.points <- rtrunc(100000, spec="gamma", a = lb.arg, b = ub.arg,
shape = shape.list[[i]], scale = scale.list[[i]])
sample.mean <- mean(sample.points)
mean.list <- append(mean.list, sample.mean)
sample.std.dev <- sd(sample.points)
std.dev.list <- append(std.dev.list, sample.std.dev)
}
Цикл for очень медленный и занимает очень много времени.
Любые лучшие решения будут оценены. Спасибо!
2 ответа
Здесь происходит несколько вещей.
Сначала вы рассчитываете весы как:
scale.list <- runif(1000, max = 100000, min = 100000)
но с тех пор min = maxвсе значения идентичны.
Во-вторых, вы не указываете lb.arg или же ub.argпоэтому я установил их на 20 и 50 произвольно.
В-третьих, профилирование этого кода с помощью Rprof показывает, что> 90% времени проводится в qtrunc(...) функция, которая вызывается rtrunc(...), Это потому, что вы генерируете 100000 точек выборки на каждой итерации, и qtrunc(...) должен сортировать это. Общее время выполнения масштабируется как O(n), где n - количество точек выборки. В моей системе использование n=1000 занимает около 7 секунд, поэтому использование n = 100000 займет 700 секунд или около 12 минут.
Мой совет: попробуйте меньше n и посмотрите, действительно ли это имеет значение. Из центральной предельной теоремы мы знаем, что распределение среднего асимптотически нормально для больших n, независимо от лежащего в основе распределения. Сомневаюсь, что увеличение n от 1000 до 100000 меняет это значительно.
Наконец, идиоматический способ сделать это в R это [использование n=1000]:
f <- function(shape,scale,lb,ub){
X <- rtrunc(1000, spec="gamma", a=lb, b=ub,shape=shape,scale=scale)
return(c(mean=mean(X),sd=sd(X)))
}
# set.seed(1) # use this for a reproducible sample
result <- mapply(f,shape.list,scale.list, lb=20,ub=50)
result.apply <- data.frame(t(result))
который создает фрейм данных с двумя столбцами: среднее и SD для каждой фигуры / масштаба. Установив начальное значение на фиксированное значение непосредственно перед запуском mapply(...)и, выполнив то же самое непосредственно перед запуском цикла for, вы можете показать, что оба они дают одинаковые результаты.
К сожалению, просто нет способа оптимизировать вашу задачу. Конечно, есть небольшие возможные оптимизации генерации случайных точек из усеченного распределения... но вот в чем дело: генерация 10^8 точек или около того из случайного распределения БУДЕТ очень медленной.
Вот несколько оптимизаций, которые я попробовал, которые немного ускоряют процесс:
генерация всех случайных чисел из равномерного распределения в [a, b] одновременно
возвращаясь к источнику определения усеченного дистрибутива, не полагаясь на "модные" пакеты (distr, distEx, truncdist)
компиляция моего кода, чтобы ускорить его
Код:
# your original code, in a function
func = function()
{
library(distr)
library(distrEx)
library(truncdist)
set.seed(42)
shape.list <- runif(1000, max = 10, min = 0.01)
scale.list <- runif(1000, max = 100000, min = 100000)
mean.list <- list()
std.dev.list <- list()
ITE.NUMBER = 10
POINTS.NUMBER = 100000
A = 0.25
B = 0.5
for (i in seq(ITE.NUMBER)) # very slow
{
sample.points <- rtrunc(POINTS.NUMBER, spec="gamma", a = A, b = B,
shape = shape.list[[i]], scale = scale.list[[i]])
sample.mean <- mean(sample.points)
mean.list <- append(mean.list, sample.mean)
sample.std.dev <- sd(sample.points)
std.dev.list <- append(std.dev.list, sample.std.dev)
}
}
# custom code
func2 = function()
{
set.seed(42)
shape.list <- runif(1000, max = 10, min = 0.01)
scale.list <- runif(1000, max = 100000, min = 100000)
mean.list <- list()
std.dev.list <- list()
ITE.NUMBER = 10
POINTS.NUMBER = 100000
A=0.25
B=0.5
#
# we generate all the random number at once, outside the loop
#
r <- runif(POINTS.NUMBER*ITE.NUMBER, min = 0, max = 1)
for (i in seq(ITE.NUMBER)) # still very slow
{
#
# back to the definition of the truncated gamma
#
sample.points <- qgamma(pgamma(A, shape = shape.list[[i]], scale = scale.list[[i]]) +
r[(1+POINTS.NUMBER*(ITE.NUMBER-1)):(POINTS.NUMBER*(ITE.NUMBER))] *
(pgamma(B, shape = shape.list[[i]], scale = scale.list[[i]]) -
pgamma(A, shape = shape.list[[i]], scale = scale.list[[i]])),
shape = shape.list[[i]], scale = scale.list[[i]])
sample.mean <- mean(sample.points)
mean.list <- append(mean.list, sample.mean)
sample.std.dev <- sd(sample.points)
std.dev.list <- append(std.dev.list, sample.std.dev)
}
}
#
# maybe a compilation would help?
#
require(compiler)
func2_compiled <- cmpfun(func2)
require(microbenchmark)
microbenchmark(func2(), func2_compiled(), func(), times=10)
Что дает следующее:
Unit: seconds
expr min lq median uq max neval
func2() 1.462768 1.465561 1.475692 1.489235 1.532693 10
func2_compiled() 1.403956 1.477983 1.487945 1.499133 1.515504 10
func() 1.457553 1.478829 1.502671 1.510276 1.513486 10
Выводы:
Как уже было сказано, улучшений практически нет: ваша задача просто очень ресурсоемкая, и с этим ничего не поделаешь.
Компиляция почти усугубила ситуацию... что и ожидается: здесь нет глупого использования плохих методов программирования (например, большие уродливые циклы)
Если вы действительно ищете улучшение скорости, вам может быть лучше с другим языком, хотя я сомневаюсь, что вы сможете добиться значительно лучших результатов..