Экспоненциальное распределение в R

Question

Я хочу смоделировать некоторые данные из распределения exp(1), но они должны быть> 0,5 . Поэтому я использовал цикл while, но, похоже, он не работает так, как мне бы хотелось. Заранее спасибо за ваши ответы!

x1<-c()

w<-rexp(1) 

while (length(x1) < 100) {

  if (w > 0.5) {

    x1<- w }

  else {

    w<-rexp(1)

  }

}

Answer 1

1) Код в вопросе имеет следующие проблемы:

• нам нужна новая случайная величина на каждой итерации, но она генерирует новые случайные величины, только если if условие ЛОЖЬ

• x1 многократно перезаписывается, а не расширяется

• хотя while может быть использован repeat кажется лучше, так как тест в конце лучше подходит, чем тест в начале

Мы можем исправить это так:

x1 <- c()
repeat {
  w <- rexp(1)
  if (w > 0.5) {
    x1 <- c(x1, w)
    if (length(x1) == 100) break
  }
}

1а) Вариант будет следующим. Обратите внимание, что if чье состояние FALSE оценивается как NULL, если нет else leg, поэтому, если условие имеет значение FALSE в строке, помеченной ##, то ничего не соединяется с x1,

x1 <- c()
repeat {
  w <- rexp(1)
  x1 <- c(x1, if (w > 0.5) w)  ##
  if (length(x1) == 100) break
}

2) Альтернативно, это генерирует 200 экспоненциальных случайных величин, сохраняя только те, которые больше, чем 0,5. Если генерируется менее 100, повторите. В конце он берет первые 100 из последней сгенерированной партии. Мы выбрали 200 достаточно большим, чтобы на большинстве прогонов требовалась только одна итерация цикла.

repeat {
  r <- rexp(200)
  r <- r[r > 0.5]
  if (length(r) >= 100) break
}
r <- head(r, 100)

Альтернатива (2) на самом деле быстрее, чем (1) или (1а), потому что она более сильно векторизована. И это несмотря на то, что оно выбрасывает больше экспоненциальных случайных величин, чем другие решения.

Answer 2

Я бы посоветовал против while (или любой другой цикл принятия / отклонения); вместо этого используйте методы из truncdist:

# Sample 1000 observations from a truncated exponential
library(truncdist);
x <- rtrunc(1000, spec = "exp", a = 0.5);

# Plot
library(ggplot2);
ggplot(data.frame(x = x), aes(x)) + geom_histogram(bins = 50) + xlim(0, 10);

Также довольно просто реализовать сэмплер, использующий выборку с обратным преобразованием, чтобы отобрать выборки из усеченного экспоненциального распределения, которое позволяет избежать отклонения выборок в цикле. Это будет более эффективный метод, чем любой метод выборки на основе принятия / отклонения, и он особенно эффективен в вашем случае, поскольку существует замкнутая форма усеченного экспоненциального cdf. Смотрите, например, этот пост для более подробной информации.