CODEWITHSUNDEEP
rregressionsimulationpoissoninverse-transform
user22130474
user22130474

Reputation:

Simulate Poisson distribution with different lambdas using inverse CDF method

I have to generate Poisson regression, for that I want to use inverse CDF method. I am following this post to generate data. But it isn't based on inverse cdf method. The post has following code

>   #sample size
> n <- 10
>   #regression coefficients
> beta0 <- 1
> beta1 <- 0.2
>   #generate covariate values
> x <- runif(n=n, min=0, max=1.5)
>   #compute mu's
> mu <- exp(beta0 + beta1 * x)
>   #generate Y-values
> y <- rpois(n=n, lambda=mu)
>   #data set
> data <- data.frame(y=y, x=x)
> data

I do not want to use the r built in function rpois. I have found this post which does utilize the method only for one lambda. Now, how can I generate Poisson distribution with inverse cdf method for a variety of lambda? N.B. This useful post provides code for variety of sample sizes, I do not want that.

Upvotes: 0

Views: 276

Answers (2)

jblood94
jblood94

Reputation: 16971

Modifying my answer in your linked post to take a vector of rates:

rpois2 <- function(n, lambda) {
  u <- runif(n)
  i <- 0L
  k <- integer(n)
  p <- exp(-lambda)
  f <- p
  idx <- 1:n
  
  while (length(idx)) {
    bln <- u[idx] < f[idx]
    k[idx[bln]] <- i
    idx <- idx[!bln]
    p[idx] <- lambda[idx]*p[idx]/(i <- i + 1L)
    f[idx] <- f[idx] + p[idx]
  }
  
  k
}

Testing:

n <- 10
beta0 <- 1
beta1 <- 0.2
x <- runif(n=n, min=0, max=1.5)
mu <- exp(beta0 + beta1 * x)
system.time({y <- replicate(1e5, rpois2(n, mu))})
#>    user  system elapsed 
#>    1.87    0.01    1.89
system.time({z <- replicate(1e5, rpois(n, mu))})
#>    user  system elapsed 
#>    0.35    0.06    0.40
data.frame(
  mean.y = rowMeans(y),
  mean.z = rowMeans(z),
  var.y = apply(y, 1, var),
  var.z = apply(z, 1, var)
)
#>     mean.y  mean.z    var.y    var.z
#> 1  2.82232 2.81855 2.826438 2.814754
#> 2  2.89834 2.91290 2.905394 2.914763
#> 3  3.41615 3.41951 3.416323 3.418016
#> 4  3.47183 3.47874 3.473361 3.460303
#> 5  3.26286 3.25541 3.261857 3.244148
#> 6  3.28565 3.28160 3.269067 3.288014
#> 7  3.49654 3.48724 3.485543 3.484132
#> 8  3.59513 3.60039 3.596586 3.604258
#> 9  3.39864 3.38836 3.378420 3.388850
#> 10 2.89590 2.89425 2.887852 2.901856

Upvotes: 1

Ben Bolker
Ben Bolker

Reputation: 226007

The answer here is elegant and fast, but if you don't care about efficiency then

my_rpois <- function(n, lambda) {
   qpois(runif(n), lambda)
}

should work to implement the inverse-CDF method. (lambda can be a vector)

Upvotes: 1

Related Questions