Оценка полученных смоделированных данных

Question

Я моделирую данные, используя метод отклонения , где функция плотности X задается как f(x)= C * e^(x) for all x in [0,1]. Я определил g(x) = 1 и C как максимум f(x), который равен 1/(e-1).

Я использовал следующий код для моделирования данных:

rejection <- function(f, C, g, rg, n) {
  naccepts <- 0
  result.sample <- rep(NA, n)

  while (naccepts < n) {
    y <- rg(1)
    u <- runif(1)

    if ( u <= f(y) / (C*g(y)) ) {
      naccepts <- naccepts + 1
      result.sample[naccepts] = y
    }
  }

  result.sample
}

f <- function(x) ifelse(x>=0 & x<=1, (exp(x)), 0)
g <- function(x) 1
rg <- runif
C <-  1/(exp(1) -1)

result <- rejection(f, C, g,rg, 1000)

Тогда Я использую histogram для сравнения смоделированных данных с curve оригинала pdf как

hist(result,freq = FALSE)
curve(f, 0, 1, add=TRUE)

Но полученный график довольно устарел! Сюжет здесь , поэтому я ищу любую помощь, чтобы уточнить, что не так в моей работе.

Answer 1

Ваше максимальное значение неверно. Для exp (x) в интервале диапазона [0 ... 1] нормализованный PDF будет

f(x) = exp(x)/(exp(1) - 1)

Таким образом, максимум f (x) равен exp (1) / (exp (1) - 1 )

ОБНОВЛЕНИЕ

Хорошо, вот код, следующий за вики-статьей о Выборка отклонения

sampleRej <- function(f, g, rg, M, n) { # rejection sampling following wiki
    naccepts <- 0
    result   <- rep(NA, n)

    while (naccepts < n) {
        y <- rg(1)
        u <- runif(1)

        if ( u <= f(y) / (M*g(y)) ) {
            naccepts <- naccepts + 1
            result[naccepts] = y
        }
    }

    result
}

g <- function(x) { # Normalized uniform distribution PDF
    1.0
}

f <- function(x) { # normalized sampling PDF
    exp(x)/(exp(1.0) - 1.0)
}

rg <- function(n) { # function to sample from g, which is uniform
    runif(n)
}


M <- exp(1.0)/(exp(1.0) - 1.0) # f(x) maximum value
q <- sampleRej(f, g, rg, M, 100000) # sample 100K points

# and plot everything together
curve(f, 0.0, 1.0)
hist(q, freq=FALSE, add=TRUE)

, и он выдает график, подобный показанному ниже, выглядит хорошо для меня

Answer 2

Это показывает всю кривую и гистограмму:

curve(f, 0, 1)
hist(result,freq = FALSE, add=TRUE)

Но, конечно, теперь гистограмма немного мала в графике ...