Pass vector of lambdas to Poisson(), or guidance on idiomatic function composition

Question

I'm trying to learn a little Julia by doing some bayesian analysis. In Peter Hoff's textbook, he describes a process of sampling from a posterior predictive distribution of a Poisson-Gamma model in which he:

1. Samples values from the gamma distribution
2. Samples values from the poisson distribution, passing a vector of lambdas

Here is what this looks like in R:

a <- 2
b <- 1

sy1 <- 217; n1 <- 111

theta1.mc <- rgamma(1000, a+sy1, b+n1)
y1.mc <- rpois(1000, theta1.mc)

In Julia, I see that distributions can't take a vector of parameters. So, I end up doing something like this:

using Distributions

a = 2
b = 1

sy1 = 217; n1 = 111

theta_mc = rand(Gamma(a+217, 1/(b+n1)), 5000)
y1_mc = map(x -> rand(Poisson(x)), theta_mc)

While I was initially put off at the distribution function not taking a vector and working Just Like R™, I like that I'm not needing to set my number of samples more than once. That said, I'm not sure I'm doing this idiomatically, either in terms of how people would work with the distributions package, or more generically how to compose functions.

Can anyone suggest a better, more idiomatic approach than my example code?

Answer 1

I would usually do something like the following, which uses list comprehensions:

a, b = 2, 1
sy1, n1 = 217, 111

theta_mc = rand(Gamma(a + sy1, 1 / (b + n1)), 1000)
y1_mc = [rand(Poisson(theta)) for theta in theta_mc]

One source of confusion may be that Poisson isn't really a function, it's a type constructor and it returns an object. So vectorization over theta doesn't really make sense, since that wouldn't construct one object, but many -- which would then require another step to call rand on each generated object.