A number of people have asked questions about problem 4 last night. Here
is a basic idea. You want to evaluate the "joint" density with different
values of y and lambda where y is poisson and lambda is gamma. This is the
joint density function, so you are not drawing random variables from any
distribution. you are simply evaluating the height of the density. To do
this, you need to create a matrix where its row and column represents y
and labmda, respectively. Then, to integrate out lambda, you want to sum
over the column. After normalizing it, this gives you the marginal
density of negative binomial.
Kosuke