Homework 7

Summary

Harvard University
Spring 2018
Instructors: Rahul Dave
**Due Date: ** Friday, March 9th, 2018 at 11:00am

Instructions:

Upload your iPython notebook containing all work to Canvas.
Structure your notebook and your work to maximize readability.

Problem 1: Gibbs Sampling On A Bivariate Normal

Let $\mathbf{X}$ be a random variable taking values in $\mathbb{R}^2$. That is, $\mathbf{X}$ is a 2-dimensional vector. Suppose that $\mathbf{X}$ is normally distributed as follows $\mathbf{X} \sim \mathcal{N} \left( \left[ \begin{array}{c} 1 \\ 2 \\ \end{array} \right], \left[ \begin{array}{ccc} 4 & 1.2 \\ 1.2 & 4 \\ \end{array} \right] \right).$ That is, the pdf of the distribution of $\mathbf{X}$ is $f_{\mathbf{X}}(\mathbf{x}) = \frac{1}{2\pi\sqrt{\vert \Sigma\vert }}\mathrm{exp}\left\{ - \frac{1}{2} (\mathbf{x} - \mu)^\top \Sigma^{-1} (\mathbf{x} - \mu)\right\}$ where $\mu = \left[ \begin{array}{c} 1 \ 2 \ \end{array} \right]$, $\Sigma = \left[ \begin{array}{ccc} 4 & 1.2 \ 1.2 & 4 \ \end{array} \right]$, and $\vert \cdot\vert $ is the matrix determinant operator.

In the following, we will denote the random variable corresponding to the first component of $\mathbf{X}$ by $X_1$ and the second component by $X_2$.

Write a Gibbs sampler for this distribution by sampling sequentially from the two conditional distributions $f_{X_1\vert X_2}, f_{X_2\vert X_1}$.
Choose a thinning parameter, burn-in factor and total number of iterations that allow you to take 10000 non-autocorrelated draws.
You must justify your choice of parameters.

Problem 2: Rubber Chickens Bawk Bawk!

In the competitive rubber chicken retail market, the success of a company is built on satisfying the exacting standards of a consumer base with refined and discriminating taste. In particular, customer product reviews are all important. But how should we judge the quality of a product based on customer reviews?

On Amazon, the first customer review statistic displayed for a product is the average rating. The following are the main product pages for two competing rubber chicken products, manufactured by Lotus World and Toysmith respectively:

Lotus World	Toysmith

Clicking on the ‘customer review’ link on the product pages takes us to a detailed break-down of the reviews. In particular, we can now see the number of times a product is rated a given rating (between 1 and 5 stars).

Lotus World	Toysmith

(The images above are also included on canvas in case you are offline, see below)

In the following, we will ask you to compare these two products using the various rating statistics. Larger versions of the images are available in the data set accompanying this notebook.

Suppose that for each product, we can model the probability of the value each new rating as the following vector: $\theta = [\theta_1, \theta_2, \theta_3, \theta_4, \theta_5]$ where $\theta_i$ is the probability that a given customer will give the product $i$ number of stars.

Part A: Inference

Suppose you are told that customer opinions are very polarized in the retail world of rubber chickens, that is, most reviews will be 5 stars or 1 stars (with little middle ground). Choose an appropriate Dirichlet prior for $\theta$. Recall that the Dirichlet pdf is given by: $f_{\Theta}(\theta) = \frac{1}{B(\alpha)} \prod_{i=1}^k \theta_i^{\alpha_i - 1}, \quad B(\alpha) = \frac{\prod_{i=1}^k\Gamma(\alpha_i)}{\Gamma\left(\sum_{i=1}^k\alpha_i\right)},$ where $\theta_i \in (0, 1)$ and $\sum_{i=1}^k \theta_i = 1$, $\alpha_i > 0 $ for $i = 1, \ldots, k$.
Write an expression for the posterior pdf, using a using a multinomial model for observed ratings. Recall that the multinomial pdf is given by: $f_{\mathbf{X}\vert \Theta}(\mathbf{x}) = \frac{n!}{x_1! \ldots x_k!} \theta_1^{x_1} \ldots \theta_k^{x_k}$ where $n$ is the total number of trials, $\theta_i$ is the probability of event $i$ and $\sum_i \theta_i = 1$, and $x_i$ is count of outcome $i$ and $\sum_i x_i = n$.

Note: The data you will need in order to define the likelihood function should be read off the image files included in the dataset.

Sample 1,000 values of $\theta$ from the posterior distribution.
Sample 1,000 values of $x$ from the posterior predictive distribution.

Part B: Ranking

Name at least two major potential problems with using only the average customer ratings to compare products.

(Hint: if product 1 has a higher average rating than product 2, can we conclude that product 1 is better liked? If product 1 and product 2 have the same average rating, can we conclude that they are equally good?)

Using the samples from your posterior distribution, determine which rubber chicken product is superior. Justify your conclusion with sample statistics.
Using the samples from your posterior predictive distribution, determine which rubber chicken product is superior. Justify your conclusion with sample statistics.
Finally, which rubber chicken product is superior?

(Note: we’re not looking for “the correct answer” here, any sound decision based on a statistically correct interpretation of your model will be fine)

Problem 3: Implementing Rat Tumors in pymc3

(it may help to see the bioassay lab to see how to structure pymc3 code, and also the examples from lecture).

Let us try to do full Bayesian inference with PyMC3 for the rat tumor example that we have solved using explicit Gibbs sampling in lab7. Remember that the goal is to estimate $\theta_i$, the probability of developing a tumor in a population of female rats that have not received treatement.

The posterior for the 70 experiments may be written thus:

$p( \{\theta_i\}, \alpha, \beta \vert Y, \{n_i\}) \propto p(\alpha, \beta) \prod_{i=1}^{70} Beta(\theta_i, \alpha, \beta) \prod_{i=1}^{70} Binom(n_i, y_i, \theta_i)$

Use uniform priors on $[0,1]$ on the alternative variables $\mu$ (the mean of the beta distribution) and $\nu$:

$\mu = \frac{\alpha}{\alpha+\beta}, \nu = (\alpha+\beta)^{-1/2}$

You may then write $\alpha$ and $\beta$ as deterministics which depend on $\mu$ and $\nu$.

Here is the data:

tumordata="""0 20 
20 
20 
20 
20 
20 
20 
19 
19 
19 
19 
18 
18 
17 
20 
20 
20 
20 
19 
19 
18 
18 
27 
25 
24 
23 
20 
20 
20 
20 
20 
20 
10 
49 
19 
46 
17 
49 
47 
20 
20 
13 
48 
50 
20 
20 
20 
20 
20 
20 
20 
48 
19 
19 
19 
22 
46 
49 
20 
20 
23 
19 
22 
20 
20 
20 
52 
46 
47 
24 
"""

Part A: Report at least the following diagostics on your samples

Autocorrelation (correlation dying by lag 20 is fine)
Parameter trace correlation after burnin
Gewecke
Gelman-Rubin
$n_{eff}$ (Number of Effective Samples)

Part B: Posterior predictive check

Recall from lab notes that in a hierarchical model there are two kinds of posterior predictions that are useful. (1) The distribution of future observations $y_i^$ given a $\theta_i$, and (2) The distribution of observations $y_j^$ drawn from a future $\theta_j$ drawn from the super-population (i.e. using the Beta on the estimated hyper parameters).

Carry out posterior predictive checks by using sample_ppc to generate posterior-predictives for all 70 experiments. This generates predictives of the first type above.
Plot histograms for these predictives with the actual value shown as a red-dot against the histogram (as in the coal disasters model in lecture 14). Is the data consistent with the predictive?

Part C: Shrinkage

Plot the posterior median of the death rate parameters $\theta_1, \theta_2, …\theta_{70}$ against the observed death rates ($y_i/n_i$)
Explain the shrinkage by comparing against a 45 degree line as done in the lab.

PART D: Experiment 71

Consider an additional experiment – experiment 71 – in which 4 out of 14 rats died.

Calculate the marginal posterior of $\theta_{71}$, the “new” experiment,
Find the $y_{71}^*$ posterior predictive for that experiment.

**HINT: ** The critical thing to notice is that the posterior including the 71st experiment factorizes:

$p(\theta_{71}, \theta_{1..70}, \alpha, \beta \vert D) \propto p(y_{71} \vert n_{71}, \theta_{71} ) p(\theta_{71} \vert \alpha, \beta) p(\theta_{1..70}, \alpha, \beta \vert D)$

Then we simply marginalize over everything to get the $\theta_{71}$ posterior:

$p(\theta_{71} \vert \theta_{1..70}, \alpha, \beta, D) = \int d\alpha \,d\beta \,d\theta_{1..70} \,p(\theta_{71}, \theta_{1..70}, \alpha, \beta \vert D)$ $= \int d\alpha \,d\beta Beta(\alpha+y_{71}, \beta + n_{71} - y_{71}) \int_{\theta_{1..70}} \,d\theta_{1..70} \,p(\theta_{1..70}, \alpha, \beta \vert D)$

The $y_{71}^*$ posterior predictive can be found in the usual way.