Goal

In our last class, we encountered a scenario we hadn’t seen before. We had not one but two parameters that we needed to estimate in our Bayesian model, which meant that our posterior distribution had two random variables in it.

When we encounter a scenario when the posterior distribution is not a known distribution (like a Gamma, a Beta, a Normal, etc.), we are not able to directly use properties of that distribution when we do our Bayesian analysis. In these situations, we need to use something called Markov Chain Monte Carlo (MCMC) to allow us to estimate the posterior distribution for our parameters of interest and hence perform inference on those parameters.

Today, we are going to explore one specific type of MCMC called Gibbs sampling. We will practice (1) deriving the full conditionals that we need, (2) coding the Gibbs sampler in R, and (3) working with the samples that we obtain from our Gibbs sampler to approximate our posterior of interest.

The Data

Today we are interested in modeling the average number of people who click on an advertisement for a new bubble tea shop in Winston-Salem. This advertisement is posted on articles for the online version of the local paper. We have information on \(Y\)= the number of people who have clicked on the ad per day for a random sample of \(n=30\) days.

Now, because \(Y\) is a count random variable, our first instinct is probably to reach for the Poisson-Gamma Bayesian model.

\[Y_1, \dots, Y_{n} \stackrel{iid}{\sim} Poisson(\theta)\]

\[\theta \sim Gamma(\alpha, \beta)\]

That’s a great starting spot! However, today we are told by our client to use something slightly different.

The Requested Bayesian Model

\[Y_1, \dots, Y_{n} \stackrel{iid}{\sim} Poisson(\theta)\]

\[\theta|\beta \sim Gamma(\alpha, \beta)\]

\[\beta \sim Exponential(c)\]

So, essentially what our client has asked us to do is to treat \(\alpha\) as a known constant. However, they do not want to set a fixed \(\beta\) hyperparameter. Instead, they want to treat \(\beta\) as a random variable.

Why would we do that?? Well, this allows us to actually model our uncertainty about \(\beta\). We are moving towards hierarchical models in this course, where we will want to be able to set priors on hyperparameters. This is a good first step towards those models!

Starting the Process

Now that we have our Bayesian model, let’s write out the posterior distribution. We have two parameters \(\theta\) and \(\beta\) that we need to estimate. This means that for Bayesian inference, we want to use the following posterior calculations:

\[P(\theta, \beta| Y_1, \dots, Y_n) = \frac{P(Y_1,\dots,Y_n| \theta,\beta)P(\theta,\beta)}{P(Y_1,\dots,Y_n)}\] Hmmm. Let’s see if we can simplify that a bit. Using properties of joint probability, we can write

\[P(\theta,\beta) = P(\theta|\beta)P(\beta)\] Also, the distribution for \(Y_i\) does not involve \(\beta\), so

\[P(Y_1,\dots,Y_n| \theta,\beta) = P(Y_1,\dots,Y_n| \theta)\]

This gives us

\[P(\theta, \beta| Y_1, \dots, Y_n) = \frac{P(Y_1,\dots,Y_n| \theta)P(\theta|\beta)P(\beta)}{P(Y_1,\dots,Y_n)}\]

We know all the parts of the numerator! However, the denominator is going to prove problematic. Even more of an issue is that fact that with two parameters, this is not going to give us a known distribution. Because of this, we will not be able to directly use properties of a known distribution to perform posterior inference.

In these settings, we need to approximate the posterior distribution through a process called Markov Chain Monte Carlo (MCMC). For today, we will use Gibbs sampling to help us approximate the posterior.

The Full Conditionals

Okay, so the goal is approximate the posterior. This means we really want to draw samples from

\[P(\theta, \beta| Y_1, \dots, Y_n) = \frac{P(Y_1,\dots,Y_n| \theta)P(\theta|\beta)P(\beta)}{P(Y_1,\dots,Y_n)}\]

We can then use these samples to do posterior inference on \(\theta\) and \(\beta\). However, we can’t directly sample from this distribution (it’s not one that is known) so we need to figure out another way to get the samples that we need.

We know that the denominator just plays the role of the normalizing constant, so

\[P(\theta, \beta| Y_1, \dots, Y_n) \propto P(Y_1,\dots,Y_n| \theta)P(\theta|\beta)P(\beta)\]

So, this is the distribution that we would like to sample from…but we can’t because there are two parameters. We can’t sample from this, but we can approximate a draw from this distribution using Gibbs Sampling.

Coding the Gibbs Sampler

Now that we have done that math that we need, let’s try coding this up! Things like Gibbs samplers tend to make more sense when we can see the algorithm running, so let’s do it.

Some last important details: We are told that in our \(n=30\) days of data, the average number of people who clicked on the ad each day was 87. The client has told us to set \(\alpha = 1125\) and \(c = .05\).

Now we are going to run just one iteration of the Gibbs sampler to start. This means we

• 1. Initialize with \(\theta = \theta_{(0)}\).
• 2. Draw \(\beta_{(1)}\) from \(\beta|\theta_{(0)}, Y_1 = y_1, \dots, Y_n =y_n\)
• 3. Draw \(\theta_{(1)}\) from \(\theta|\beta_{(1)}, Y_1 = y_1, \dots, Y_n =y_n\)

This is the core of the code that we need! Now that we can do one iteration, let’s perform our Gibbs sampler for \(S = 10000\) iterations.

When you finish Question 9, you should have an 10000 by 2 data frame filled with your posterior draws for \(\theta\) and \(\beta\). Now what??

Assessing Convergence

We now have the results of our Gibbs sampler - fantastic!! However, there is one thing that we need to do before we can use our results. We need to pause and assess convergence. This means we need to check to make sure that the results of our draws look fairly consistent from sample to sample.

We don’t want the results to be identical (we need to see the variability from sample to sample in order to estimate the distribution of the parameter). However, wildly different values for each sample can mean that we need to run the process longer.

As we discussed in class, one method that we use to assess convergence is a trace plot.

It is important to check for convergence before we proceed with any posterior inference. If you do not check, and the sampler has not converged, your draws from your full conditionals will not approximate draws from your posterior of interest.

There are two other things that we are going to explore before we start using our posterior samples.

Using the Samples

At this point, we have explored convergence under a variety of settings. Let’s return to our original Gibbs sampler (when we initialized \(\theta_{(0)}\)) and do some posterior analysis!

We have posterior samples for \(\beta\) and \(\theta\), and we could do inference on both. \(\theta\) is really the parameter of interest, but let’s explore both for now.

suppressMessages(library(HDInterval))
HDInterval::hdi( GibbsSamples$beta, credMass=.95)