如何修改从pdf生成随机数的功能更快运行

How to modify the function of generating random numbers from a pdf to run faster

下面是从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)
    }

但是,该函数花了很长时间来模拟随机变量,有没有办法将此代码修改为更快运行?

您的代码中存在一些错误,导致您的实施需要很长时间。

密度函数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))
}

请注意,您需要将 lgamma 加 1 作为 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)
    }

希望对您有所帮助!