Экспоненциальное распределение в 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 может использоваться * Кажется 1018 *лучше, поскольку тест в конце подходит лучше, чем тест в начале

Мы можем исправить это следующим образом:

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

1a) Вариант будет следующим.Обратите внимание, что if с условием FALSE оценивается как NULL, если отсутствует нога else, поэтому, если условие имеет значение 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.См. Например этот пост для более подробной информации.