лучший способ создать функцию, которая возвращает гамма-распределение?

Поэтому я создал функцию, которая генерирует кадр данных для гамма-распределения.

В настоящее время у меня есть.

sample_gamma <- function(alpha,beta,n,iter) {
gamma.df <- as.data.frame(matrix(nrow = iter, ncol = 3))
colnames(gamma.df) <- c("iteration","mean","standard dev")
gamma.df$iteration <- c(1:iter)
  for (i in 1:iter) {
    gamma.dist <- rgamma(n,shape = alpha, rate = beta, scale = 1/beta)
    gamma.df[i,2] <- mean(gamma.dist)
    gamma.df[i,3] <- sd(gamma.dist)
  }
print(gamma.df)
}

Функция делает все, что мне нужно, но мне было интересно, есть ли какие-либо альтернативные или более чистые способы сделать это.


r function statistics distribution gamma-distribution
person Gian Batayola    schedule 21.11.2019    source источник
comment
Должен быть класс, который изучает это, как я нашел stackoverflow.com/q/58964331, stackoverflow.com/q/58955845 и stackoverflow.com/q/58955845, которые генерируют функцию с именем sample_gamma, включающую одни и те же формальные аргументы и близкие к одному и тому же результату.   -  person r2evans    schedule 21.11.2019

Ответы (2)


arrow_upward
3
arrow_downward

Я бы создал функцию, которая возвращает mean и sd за одну итерацию.

sample_gamma <- function(alpha,beta,n) {
    dist <- rgamma(n,shape = alpha, rate = beta)
    c(mean = mean(dist), sd = sd(dist))
}

а затем повторите это, используя replicate

t(replicate(5, sample_gamma(2, 3, 4)))

#          mean        sd
#[1,] 0.5990206 0.2404226
#[2,] 0.6108976 0.3083426
#[3,] 1.0616542 0.4602403
#[4,] 0.3415355 0.1543885
#[5,] 1.0558066 0.9659599
person Ronak Shah    schedule 21.11.2019
comment
Лучший ответ на все остальные (иначе идентичные) вопросы. - person r2evans; 21.11.2019

arrow_upward
2
arrow_downward

Хотя я думаю, что ответ Ронака Шаха прост и относительно идиоматичен (с точки зрения R), вот ответ, который немного более эффективен при масштабировании до высоких значений iter (поскольку он выполняет только один случайный выбор):

sample_gamma <- function(alpha, beta, n, iter) {
  mtx <- matrix(rgamma(n*iter, shape=alpha, rate=beta), nrow=n, ncol=iter)
  t(apply(mtx, 2, function(a) c(mean=mean(a), sd=sd(a))))
}
sample_gamma(2, 3, 4, 5)
#           mean         sd
# [1,] 0.6486220 0.22900833
# [2,] 0.8551055 0.07874287
# [3,] 0.7854750 0.72694260
# [4,] 0.7045878 0.24834502
# [5,] 1.1783301 0.25210538

Сравнительный анализ:

microbenchmark::microbenchmark(
  RS=t(replicate(5, sample_gamma_RS(2,3,4))),
  r2=sample_gamma_r2(2,3,4,5)
)
# Unit: microseconds
#  expr   min     lq    mean median    uq    max neval
#    RS 413.7 493.70 757.884 743.80 946.1 1611.6   100
#    r2 405.2 461.15 681.630 706.35 898.6 1348.2   100

microbenchmark::microbenchmark(
  RS=t(replicate(500, sample_gamma_RS(2,3,4))),
  r2=sample_gamma_r2(2,3,4,500)
)
# Unit: milliseconds
#  expr    min       lq     mean   median       uq      max neval
#    RS 31.271 40.58735 56.44298 57.85735 65.08605  95.1866   100
#    r2 29.110 38.81230 53.99426 57.45820 61.35720 100.5820   100

microbenchmark::microbenchmark(
  RS=t(replicate(500, sample_gamma_RS(2,3,400))),
  r2=sample_gamma_r2(2,3,400,500)
)
# Unit: milliseconds
#  expr     min       lq     mean   median       uq      max neval
#    RS 60.6782 101.3112 121.3533 116.7464 140.8845 227.1904   100
#    r2 66.3892  81.0329 106.9920  98.7170 126.7742 198.3947   100

Признаюсь, я думал, что разница в производительности будет более существенной.

person r2evans    schedule 21.11.2019
comment
Тоже хороший подход. Если вы хотите рассчитать среднее значение по столбцам и sd, вы можете использовать base::colMeans и matrixStats::colSds, что должно улучшить производительность. - person Ronak Shah; 21.11.2019
comment
Правда, эти функции могут выжать из него немного больше скорости. Почему-то я не уверен, что класс, изучающий этот sample_gamma, беспокоится о производительности больших выборок, но это было занимательное упражнение. :-) - person r2evans; 21.11.2019