13.8 Compound Distribution (Optional)

An example of where simulation is useful is shown below. Suppose you have one variable, \(X\) that is assumed to have a Poisson distribution with lambda = 3 and another random variable \(Y\) that is assumed to have a Gamma distribution with shape = 3 and rate = 0.5.

A practical story where this scenario might be useful is that we want to explore the annual sum of costs of individuals with car accidents. The number of annual car accidents is Poisson and the cost per car accident is Gamma.

The code below contains a for-loop as discussed in Section 7.4 and data frames as discussed in Section 10.

The example assumes that there are 6 people to consider. The data frame is initialized in the first block of code where one Poisson is simulated for the number of accidents, Gamma values are generated for the number of accidents, and then summed over the total. The number of accidents and the sum is written to a data frame. A second block of code contains a for loop that does the same thing, but adds the results to the already created data frame.

The interim values x, y, z are printed so that you can verify the process and the final results.

  x <- rpois(n = 1, lambda = 3)
  y <- rgamma(n = x, shape = 3, rate = 0.5) 
  z <- sum(y)
  print(x)
  print(y)
  print(z)
  dat <- data.frame(x,z)
  dat

for (count in 1:5){
  x <- rpois(n = 1, lambda = 3)
  y <- rgamma(n = x, shape = 3, rate = 0.5) 
  z <- sum(y)
  print(x)
  print(y)
  print(z)
  dat[nrow(dat) + 1,] <- list(x, z)
       }

}

If the simulation is completed with more iterations, we can plot the new simulated distribution, along with some summary statistics.

  x <- rpois(n = 1, lambda = 3)
  y <- rgamma(n = x, shape = 3, rate = 0.5) 
  z <- sum(y)
  dat1 <- data.frame(x,z)
  

for (count in 1:999){
  x <- rpois(n = 1, lambda = 3)
  y <- rgamma(n = x, shape = 3, rate = 0.5) 
  z <- sum(y)
  dat1[nrow(dat1) + 1,] <- list(x, z)
}
  
  hist(dat1$z)
  describe(dat1)