Как изменить функцию генерации случайных чисел из PDF, чтобы работать быстрее

Question

Ниже приведен код генерации случайных чисел из pdf:

Я модифицирую код из функции rcmp (пакет COMPoissonReg).

  dcomp <- function(y,mu,v,z=NULL, max=100)
    {
      if (is.null(z)){
      z=sum(sapply(  0:100, function(j) (( ((mu)^j) / (factorial(j)) )^v)  ))
      }
      log.ff <- v*y*log(mu)-v*lgamma(y) - log(z)
      return(exp(log.ff))
    }

    rcomp <- function(n, mu, v, max=100)
    {
      if (length(mu) == 1) {
        mu <- rep(mu, n)
      }
      if (length(v) == 1) {
        v <- rep(v, n)
      }
      u <- runif(n)
      y <- numeric(n)
      z=sum(sapply(  0:100, function(j) (( ((mu)^j) / (factorial(j)) )^v)  ))
      for (i in 1:n) {
        px <- dcomp(y[i], mu[i], v[i],z=z[i], max = max)
        while (px < u[i]) {
          y[i] <- y[i] + 1
          px <- px + dcomp(y[i], mu[i], v[i], z=z[i],max = max)
        }
      }
      return(y)
    }

Однако функции имитировали случайные переменные очень долго, есть ли способ изменить этот код, чтобы он работал быстрее?

Answer 1

В вашем коде есть ряд ошибок, из-за которых ваша реализация занимает много времени.

Функция плотности dcomp должна быть изменена следующим образом

dcomp <- function(y,mu,v,z=NULL, max=100) {
  if (is.null(z)){
    z=sum(sapply(  0:100, function(j) (( ((mu)^j) / (factorial(j)) )^v)  ))
  }
  log.ff <- v*y*log(mu)-v*lgamma(y+1) - log(z)
  return(exp(log.ff))
}

Примечаниечто вам нужно добавить 1 к lgamma как gamma (x + 1) = factorial (x).

В функции rcomp, где вы генерируете случайные величины, у вас есть проблема в

• sum(sapply( 0:100, function(j) (( ((mu)^j) / (factorial(j)) )^v) ))

  sum в этой строке сворачивает векторизацию.Вам необходимо обновить его, чтобы получить правильный вектор с индивидуализированными значениями z.Я удалил этот предварительный расчет в приведенном ниже коде и просто выполняю вычисления внутри dcomp, но выполняю предварительный расчет с определенно экономящим время.

Обновленная функция rcomp выглядит следующим образом

rcomp <- function(n, mu, v, max=100)
{
  if (length(mu) == 1) {
        mu <- rep(mu, n)
      }
      if (length(v) == 1) {
        v <- rep(v, n)
      }
      u <- runif(n)
      y <- rep(0, n)  # Have changed this line to force zeros as starting points
      # z=sum(sapply(  0:100, function(j) (( ((mu)^j) / (factorial(j)) )^v)  ))
      for (i in 1:n) {
        px <- dcomp(y[i], mu[i], v[i], max = max) # Not using z
        while (px < u[i]) {
          y[i] <- y[i] + 1
          px <- px + dcomp(y[i], mu[i], v[i],max = max)  # Also not using z
        }
      }
      return(y)
    }

Надеюсь, это поможет!