3
\$\begingroup\$

I have an exercise that ask to plot the prior and posterior distribution in the Poisson-Gamma model. I did it (I think it's correct) inspired in the answer of this question https://stats.stackexchange.com/questions/70661/how-does-the-beta-prior-affect-the-posterior-under-a-binomial-likelihood

Is there anything that can improve the code?

colors = c("red","blue","green","orange","purple")

n = 10

lambda = .2

x = rpois(n,lambda)
grid = seq(0,7,.01)


alpha = c(.5,5,1,2,2)
beta = c(.5,1,3,2,5)

plot(grid,grid,type="n",xlim=c(0,5),ylim=c(0,4),xlab="",ylab="Prior Density",
     main="Prior Distributions", las=1)
for(i in 1:length(alpha)){
  prior = dgamma(grid,shape=alpha[i],rate=1/beta[i]) 
  lines(grid,prior,col=colors[i],lwd=2)
}

legend("topleft", legend=c("Gamma(0.5,0.5)", "Gamma(5,1)", "Gamma(1,3)", "Gamma(2,2)", "Gamma(2,5)"),
       lwd=rep(2,5), col=colors, bty="n", ncol=3)

for(i in 1:length(alpha)){
  dev.new()
  plot(grid,grid,type="n",xlim=c(0,5),ylim=c(0,10),xlab="",ylab="Density",xaxs="i",yaxs="i",
       main="Prior and Posterior Distribution")

  alpha.star = alpha[i] + sum(x)
  beta.star = beta[i] + n
  prior = dgamma(grid,shape=alpha[i],rate=1/beta[i])
  post = dgamma(grid,shape=alpha.star,rate=beta.star)

  lines(grid,post,lwd=2)
  lines(grid,prior,col=colors[i],lwd=2)
  legend("topright",c("Prior","Posterior"),col=c(colors[i],"black"),lwd=2)

}
\$\endgroup\$
  • \$\begingroup\$ Could I ask why you use 1/beta[i] in the first set of priors and beta[i] in the second set of priors.The first thing I'd do if refactoring your code would be to write a function that does the code that you've written inside the for-loops; then I'd split the data-generation from the data-plotting or plot-appending logic into different functions. If I was rethinking how you present your plots, I'd consider plotting several prior/posterior figures side-by-side, eg, using par(mfrow = c(1, 5)) and add a meaningful x-label \$\endgroup\$ – Russ Hyde Jan 25 at 10:26
  • \$\begingroup\$ I'd also expand the range of the x-axis in your first plot: < 20% of the density is covered for all but the Gamma(0.5, 0.5) distibution \$\endgroup\$ – Russ Hyde Jan 25 at 10:32
  • 1
    \$\begingroup\$ I have another query: your choices of alpha and beta jump around quite a bit. The mean of the gamma-dist defined by your alpha/beta pairs varies between 0.3 and 5, and the variance from 0.08 to 5. It might be more illustrative to pick three prior means and two prior variances (say); and then choose alpha/beta pairs that are consistent with the 6 possible combinations of mean/variance \$\endgroup\$ – Russ Hyde Jan 25 at 11:10
  • \$\begingroup\$ @RussHyde ooh, it's a typo, should be 1/beta[i]. With this prior/posterior figures side-by-side you mean one prior in one plane and the posterior in a different plane but next to the prior? \$\endgroup\$ – user178403 Jan 25 at 19:21
  • \$\begingroup\$ @RussHyde I've change to xlim=c(0,2) now looks complete \$\endgroup\$ – user178403 Jan 25 at 19:22
3
\$\begingroup\$

There was a couple of errors in the code as originally posted. The code was OK as far as base R goes.

First thing I did was run styler and then lintr on your code; these two things help clean up the coding style in your scripts.

That does things like this:

colors = c("red","blue","green","orange","purple")

# changed to (spaces / idiomatic assignment):
colors <- c("red", "blue", "green", "orange", "purple")

Then I changed your 1:length(alphas) to seq_along(alphas). The latter is a bit safer since the former can fail with empty input.

Then I replaced your 5-separate prior/posterior plots with a single plot that contains 5 panels. This makes it easier to compare the appropriateness of the different priors. To do this, I removed your dev.new()s, added a call to par(mfrow = c(number_of_rows, number_of_columns)) and obviously, tidied this up afterwards (returning to a 1x1 grid)

par(mfrow = c(2, ceiling(length(alpha) / 2)))

for (i in seq_along(alpha)) {
    # removed dev.new()
   ... plotting code ...
   )
 }

par(mfrow = c(1, 1))

Then I cleaned up your experimental data and your prior-parameters / plotting parameters; they were all winding around each other. I also renamed your alpha / beta vectors - in R, these correspond to the shape and rate parameters that are passed into dgamma:

# ---- experimental data

num_observations <- 10

lambda <- .2

x <- rpois(num_observations, lambda)

# ---- prior parameters

# assumed 'beta' was a rate parameter
# - this, since there was confusion in the parameterisation of dgamma():
#   - early section used rate = 1 / beta[i];
#   - later section used rate = beta[i]; and
#   - definition of beta_star = beta[i] + n; implied beta was definitely a rate

shape <- c(.5, 5, 1, 2, 2)
rate <- c(.5, 1, 3, 2, 5)

# ---- plotting parameters

colors <- c("red", "blue", "green", "orange", "purple")

# ---- search parameters

grid <- seq(0, 2, .01)

Then I made a function to do your prior-comparison stuff (the first set of plots). Any parameters that were to be passed through to plot were passed in using the ... argument.

# ---- comparison of the prior distributions

plot_priors <- function(grid, shapes, rates, colors,
                        legend_text, lwd = 2, ...) {
  plot(grid, grid, type = "n", ...)

  for (i in seq_along(shape)) {
    prior <- dgamma(grid, shape = shape[i], rate = rate[i])
    lines(grid, prior, col = colors[i], lwd = lwd)
  }

  legend(
    "topleft",
    legend = legend_text, lwd = lwd, col = colors, bty = "n", ncol = 2
  )
}

This can be called like:

plot_priors(
  grid, shape, rate, colors,
  legend_text = paste0("Gamma(", c("0.5,0.5", "5,1", "1,3", "2,2", "2,5"), ")"),
  xlim = c(0, 1), ylim = c(0, 4), xlab = "", ylab = "Prior Density",
  main = "Prior Distributions", las = 1
)

It's useful to split your computations away from your plotting code - so I extracted the code you used to compute the posterior params:

compute_posterior_parameters <- function(observations,
                                         prior_shape,
                                         prior_rate) {
  list(
    shape = prior_shape + sum(observations),
    rate = prior_rate + length(observations)
  )
}

Then I pulled the plotting code for your prior/posterior comparisons into a function (similarly to the above)

plot_prior_post_comparison <- function(
                                       observations,
                                       grid, shapes, rates, colors,
                                       lwd = 2,
                                       ...) {
  # make a grid for plotting
  par(mfrow = c(2, ceiling(length(shapes) / 2)))

  for (i in seq_along(shapes)) {
    # details of the prior and post distributions
    posterior_params <- compute_posterior_parameters(
      observations,
      prior_shape = shapes[i], prior_rate = rates[i]
    )
    prior <- dgamma(
      grid,
      shape = shapes[i],
      rate = rates[i]
    )
    post <- dgamma(
      grid,
      shape = posterior_params$shape,
      rate = posterior_params$rate
    )

    # plotting code
    plot(grid, grid, type = "n", ...)
    lines(grid, post, lwd = lwd)
    lines(grid, prior, col = colors[i], lwd = lwd)
    legend("topright",
      c("Prior", "Posterior"),
      col = c(colors[i], "black"), lwd = lwd
    )
  }

  # revert the plotting grid back to 1x1
  par(mfrow = c(1, 1))
}

note that the par calls, which change the plotting grid are all nested inside the function, so any subsequent plots should be unaffected.

Then I called that function:

# ---- prior/posterior comparison

plot_prior_post_comparison(
  observations = x,
  grid = grid, shapes = shape, rates = rate, colors = colors,
  xlim = c(0, 1), ylim = c(0, 10), xlab = "", ylab = "Density",
  xaxs = "i", yaxs = "i",
  main = "Prior and Posterior Distribution"
)


Then I put all the functions at the start and all the calls at the end of a script:

The full code:

# ---- comparison of the prior distributions

plot_priors <- function(grid, shapes, rates, colors,
                        legend_text, lwd = 2, ...) {
  plot(grid, grid, type = "n", ...)

  for (i in seq_along(shape)) {
    prior <- dgamma(grid, shape = shape[i], rate = rate[i])
    lines(grid, prior, col = colors[i], lwd = lwd)
  }

  legend(
    "topleft",
    legend = legend_text, lwd = lwd, col = colors, bty = "n", ncol = 2
  )
}

# ---- prior:posterior analysis

compute_posterior_parameters <- function(observations,
                                         prior_shape,
                                         prior_rate) {
  list(
    shape = prior_shape + sum(observations),
    rate = prior_rate + length(observations)
  )
}

plot_prior_post_comparison <- function(
                                       observations,
                                       grid, shapes, rates, colors,
                                       lwd = 2,
                                       ...) {
  # make a grid for plotting
  par(mfrow = c(2, ceiling(length(shapes) / 2)))

  for (i in seq_along(shapes)) {
    # details of the prior and post distributions
    posterior_params <- compute_posterior_parameters(
      observations,
      prior_shape = shapes[i], prior_rate = rates[i]
    )
    prior <- dgamma(
      grid,
      shape = shapes[i],
      rate = rates[i]
    )
    post <- dgamma(
      grid,
      shape = posterior_params$shape,
      rate = posterior_params$rate
    )

    # plotting code
    plot(grid, grid, type = "n", ...)
    lines(grid, post, lwd = lwd)
    lines(grid, prior, col = colors[i], lwd = lwd)
    legend("topright",
      c("Prior", "Posterior"),
      col = c(colors[i], "black"), lwd = lwd
    )
  }

  # revert the plotting grid back to 1x1
  par(mfrow = c(1, 1))
}


# ----

# ---- experimental data

num_observations <- 10

lambda <- .2

x <- rpois(num_observations, lambda)

# ---- prior parameters

# assumed 'beta' was a rate parameter
# - this, since there was confusion in the parameterisation of dgamma():
#   - early section used rate = 1 / beta[i];
#   - later section used rate = beta[i]; and
#   - definition of beta_star = beta[i] + n; implied beta was definitely a rate

shape <- c(.5, 5, 1, 2, 2)
rate <- c(.5, 1, 3, 2, 5)

# ---- plotting parameters

colors <- c("red", "blue", "green", "orange", "purple")

# ---- search parameters

grid <- seq(0, 2, .01)

# ---- comparison of priors

plot_priors(
  grid, shape, rate, colors,
  legend_text = paste0("Gamma(", c("0.5,0.5", "5,1", "1,3", "2,2", "2,5"), ")"),
  xlim = c(0, 1), ylim = c(0, 4), xlab = "", ylab = "Prior Density",
  main = "Prior Distributions", las = 1
)

# ---- prior/posterior comparison

plot_prior_post_comparison(
  observations = x,
  grid = grid, shapes = shape, rates = rate, colors = colors,
  xlim = c(0, 1), ylim = c(0, 10), xlab = "", ylab = "Density",
  xaxs = "i", yaxs = "i",
  main = "Prior and Posterior Distribution"
)

I don't thnk its perfect - the plotting and calculation steps are still pretty tied together - but it should be more easy to add extra pairs of rate/shape values into your prior distributions, for example.

One place the code coould be further improved is by passing in a data.frame where each row contains the shape/rate/colour values and any annotations for a given prior. The functions contain too many arguments at the moment, and this would fix that (and guarantee that there is a shape for each rate).

\$\endgroup\$
  • \$\begingroup\$ Thank you so much Russ Hyde. I've tested the full code and the output looks really nice :) . I'll read in detail your full answer later to understand the code. \$\endgroup\$ – user178403 Jan 25 at 20:46
  • \$\begingroup\$ With this: the plotting and calculation steps are still pretty tied together, you meant the code could be improved? Or you meant the output of the code could be improved ? Or both? \$\endgroup\$ – user178403 Jan 26 at 18:42
  • \$\begingroup\$ More the code than the output. It's important to think about what bits of the code you would need to change if you wanted to do X. X could be any number of things: i) you might want to sample from a different class of prior; ii) you might want to add an extra example plot or change the parameters for your current examples; iii) you might want to plot using ggplot2 instead of base etc. At present, to do any of these things you'd have to change a few different places in the code, so there's still improvements that could be made \$\endgroup\$ – Russ Hyde Jan 27 at 11:57
  • \$\begingroup\$ I see. For example, What new parameters would be a good idea to use? \$\endgroup\$ – user178403 Jan 27 at 19:11
  • \$\begingroup\$ Haha, that's for you to decide. I was just giving an example of how the code might evolve... \$\endgroup\$ – Russ Hyde Jan 28 at 10:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.