---
title: "Lab8 - sampling"
author: "Lukasz Kozlowski (based on Michał Ciach meterial)"
date: "December 06, 2018"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Monte Carlo approach to Chocolate

Today we will use more advanced techniques to analyze the chocolate dataset. Recall that in the previous exercises (lab6), we've assumed that the prior probability of being tasty was $Beta(a, b)$. We've parameterized this distribution with $n=a+b$ and $m = a/(a+b)$, i.e. the number of chocolate bars in our prior knowledge and the assumed average tastiness. In this parametrization, the prior was $Beta(mn, (1-m)n)$. The parameter $n$ represented the strength of our prior knowledge - in other words, it was equal to the number of chocolate bars that we would need to inspect so that we would equally believe in the data and our prior knowledge. 

We assumed that $90\%$ bars are tasty, and one would need 20 bars to change his/her beliefs. This lead to $m=0.9$ and $n=20$. In reality, we are not exactly sure if this should be $n=20$ - maybe it's $21$ or $18$? In fact, it's reasonable to assume that $m$ and $n$ are random as well. This leads us to the notion of *hierarchical models*, which have multiple layers of randomness.

As before, let $X$ be the number of tasty chocolate bars, and assume that $X \sim Bin(n, \theta)$ and $\theta \sim Beta(nm, n(m-1))$. This time, however, we'll assume that $n$ and $m$ are random. It's reasonable to assume that $n$ is from Gamma distribution (which is a bit similar to Poisson distribution, but continuous), and $m$ from Beta. The full specification of our model is as follows:

\begin{align}
X & \sim Bin(n, \theta), \\
\theta & \sim Beta(nm, n(1-m)), \\
n & \sim Gamma(s, r), \\
m & \sim Beta(\mu, \nu).
\end{align}

The deterministic variables $s, r, \mu$ and $\nu$ are called *hyperparameters*, because they influence the *parameter* $\theta$. The values of the hyperparameters represent our prior knowledge. Our goal is to infer the **posterior** distributions of $\theta, n$ and $m$. This is, however, usually impossible. Instead, we'll use the Monte Carlo approach to simulate the random parameters from the posterior distribution. 

There are many methods to sample from the posterior distributions, the most important being the Metropolis-Hastings algorithm and the Gibbs sampler. The latter often has a better performance, but also often it's very difficult to use. Therefore, we'll start with Metropolis-Hastings (MH). 

The theory between the MH algorithm has been covered during lecture, and I assume you more or less know the idea behind it. This algorithm, in it's basic form, is very simple, and it's a very good exercise to implement it on your own. However, we won't do that here; Instead, we'll focus on learning STAN, a programming language for hierarchical bayesian modelling. STAN has interfaces for many popular languages, including R and Python. Here, we'll focus on the R interface; You can use Python if you prefer, but you'll have to install it by yourself. Apart from installation, using STAN is exactly the same in all languages.


## STAN

Adapted from [this tutorial](http://nbviewer.jupyter.org/github/QuantEcon/QuantEcon.notebooks/blob/master/IntroToStan_basics_workflow.ipynb). Full documentation of STAN language is [here](http://mc-stan.org/users/documentation/). 

First, install the `rstan` package. STAN is written in C++, so the installation needs to take some time to compile it.

```{r, eval=FALSE}
install.packages('rstan')
```


Load the installed package. We'll also load some additional packages and enable multithreading.
 
```{r}
library(ggplot2); library(rstan); library(reshape2)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
```

A Stan model is comprised of code blocks. Each block is a place for a certain task. The bold blocks below must be present in all Stan programs (even if they contain no arguments):

  1. functions, where we define functions to be used in the blocks below. This is where we will write out a random number generator that gives us draws from our assumed model.
  2. data, declares the data to be used for the model
  3. transformed data, makes transformations of the data passed in above
  4. parameters, defines the unknowns to be estimated, including any restrictions on their values.
  5. transformed parameters, often it is preferable to work with transformations of the parameters and data declared above; in this case we define them here.
  6. model, where the full probability model is defined.
  7. generated quantities, generates a range of outputs from the model (posterior predictions, forecasts, values of loss functions, etc.).

First, the model, composed of the above blocks, needs to be written as a string. Below is an example of a model written in STAN language; this is the model from the previous tutorial. 

```{r}
chocolate.model <- "
data {
  // In this section, we define the data that must be passed to Stan (from whichever environment you are using)

  int N; // number of tested chocolate bars
  int X; // number of tasty bars
  real n;
  real m; // hyperparameters n and m
}
parameters {
  // Define the parameters that we will estimate, as well as any restrictions on the parameter values (e.g. standard deviations can't be negative...)
  real<lower=0, upper=1> theta;
}
transformed parameters {
  vector[2] beta_pars;  // vector of parameters
  beta_pars[1] = n*m;
  beta_pars[2] = n*(1-m);
}
model {
  // This is where we write out the probability model, in very similar form to how we would using paper and pen

  // Define the priors
  theta ~ beta(beta_pars[1], beta_pars[2]); // a= n*m, b=n*(1-m)
  
  // The likelihood
  X ~ binomial(N, theta);
}
generated quantities {
  // All the parameters are simulated by default, so no need to write them here. We will also generate posterior predictive draws (yhat) for each data point. These will be elaborated on below. 
  
  int y_sim;
  y_sim = binomial_rng(N, theta);
}"
```

Now, load the data and set the prior hyperparameters. 

```{r}
library(readr)
chocolate <- read_csv("chocolate.csv")
chocolate$Good <- ifelse(chocolate$Rating > 2, "Yes", "No")
prior.n <- 20
prior.m <- 0.9
#chocolate$Good 
```

Specify the data list that we will pass to Stan. This gives Stan everything declared in the `data{}` block. 

```{r}
data_list <- list(N = nrow(chocolate), X = sum(chocolate$Good=="Yes"), n = prior.n, m = prior.m)
```


Call Stan. You'll need to specify either `model_code` (like the ones we defined above), a `file` (.stan file, containing code), or a fitted Stan object (`fit`). You should also pass Stan a data list, number of cores to estimate on, the number of Markov chains to run (4 by default) and number of iterations (2000 by default). 

We'll run just one chain with 12 000 iterations. Here, an iteration is a step of a Markov chain. Each step realizes sampling from the posterior distribution. We will discard the first 2000 iterations. This is called *burn-in* or *warm-up*, and is used to let the chain converge to the stationary distribution. 

The first time you run the model, it needs to compile the code; this may take some time. 

```{r}
chocolate.fit <- stan(model_code = chocolate.model, data = data_list, warmup=2000, chains = 1, iter = 12000)
```

After STAN has run the simulation, we can extract the results.

```{r}
thetas <- extract(chocolate.fit, "theta", permuted=F)
posteriors <- extract(chocolate.fit, "y_sim", permuted=F)
```

Finally we can visualize the posterior distribution of $\theta$, and perform the *posterior predictive check* - i.e. visualize the distribution of the data under the posterior distribution of the parameter. 

First, compare the simulation results with theoretical predictions.

```{r}
qplot(c(thetas), binwidth=0.001, xlim=c(0.0,1))
```

```{r}
qplot(seq(0, 1, by=0.001), dbeta(seq(0, 1, by=0.001), 
                                 shape1=prior.n*prior.m + sum(chocolate$Good=="Yes"), 
                                 shape2=prior.n*(1-prior.m) + nrow(chocolate) - sum(chocolate$Good=="Yes")),
      geom="line")
```


Next, compare the posterior data distributions. 

```{r}
qplot(c(posteriors), xlim=c(1600, 1800), binwidth=1, main="Posterior distribution from Stan Monte Carlo simulation")
```

```{r}
pred.post.X <- rbinom(10000, nrow(chocolate), rbeta(10000, 
                                                    shape1=prior.n*prior.m + sum(chocolate$Good=="Yes"), 
                                                    shape2=prior.n*(1-prior.m) + nrow(chocolate) - sum(chocolate$Good=="Yes")))
qplot(pred.post.X, xlim=c(1600, 1800), binwidth=1, main="Posterior distribution from our own simulations")
```


Both plots seem exactly the same! This means that our method works, and our STAN model was properly specified. Let's now use STAN to simulate posterior distributions in our simple hierarchical model.

# Simple hierarchical model

Recall that, at the beginning of this tutorial, we have discussed the following hierarchical model:

\begin{align}
X & \sim Bin(n, \theta), \\
\theta & \sim Beta(nm, n(1-m)), \\
n & \sim Gamma(s, r), \\
m & \sim Beta(\mu, \nu).
\end{align}

The only difference from the previous model is that now we want to estimate the posterior distributions of all three parameters: $\theta, n$ and $m$. As before, $n$ measures the influence of $m$ on $\theta$, but this time these variables are random. Deriving analytical equations for the distributions of parameters was simple in the previous model, but would be very difficult now. This is where MCMC methods really come in handy.

Our model is parameterized by hyperparameters $s, r, \mu$ and $\nu$. Those hyperparameters, together with the total number of tasty chocolate bars and number of bars tested, will be the arguments of the STAN model. 

To obtain easily interpretable parameters for $Gamma(s, r)$, note that if $Y \sim Gamma(s, r)$, then

\begin{align}
\mathbb{E}X &= s/r \\
Var X & = s/r^2
\end{align}

It follows that

\begin{align}
s & = \mathbb{E}^2 X / Var X \\
r & = \mathbb{E} X / Var X.
\end{align}

Now, recall that $n$ represents the strength of our prior knowledge (the number of chocolate bars we'd need to check to alter our beliefs). So, if we believe that our beliefs are worth about $20 \pm 10$ chocolate bars, it's reasonable to parametrize the prior by specifying

\begin{align}
s & = 10^2 / 10^2 = 1 \\
r & = 10 / 10^2 = 0.1.
\end{align}

**Your task:**

* Plot the density of the prior Gamma for $s=4, r=0.2$ on a line plot.  

    * Use `dgamma` function and either `qplot(x, y, geom="line")` or `plot(x, y, type="l")`.
    
    * An example interval over which to plot is $[0, 100]$.
    
```{r, include=F}
qplot(seq(0, 100, by=0.2), dgamma(seq(0, 100, by=0.2), shape = 4, rate=0.2), geom="line", ylab="f(x)", xlab="x")
```

* Plot the density of the prior Beta for $\mu =4$ and $\nu=2$.

```{r, include=F}
qplot(seq(0, 1, by=0.001), dbeta(seq(0, 1, by=0.001), 4, 2), geom="line", ylab="f(x)", xlab="x")
```


* Implement the above hierarchical model in STAN.

    * Use the previous model
    * Remember about the restrictions for parameters (e.g. $n \geq 0$)

* Run the analysis for $s=4, r=0.2, \mu = 4$ and $\nu=2$. You may as well pick your own hyperparameters. Use one chain and 4 000 iterations, including 2 000 warmup iterations.



```{r, include=F}
chocolate.hierarchical.model <- "
data {
  // In this section, we define the data that must be passed to Stan (from whichever environment you are using)

  int N; // number of tested chocolate bars
  int X; // number of tasty bars
  real s;
  real r; // hyperparameter of n 
  real mu;
  real nu; // hyperparameters of m
}
parameters {
  // Define the parameters that we will estimate, as well as any restrictions on the parameter values (e.g. standard deviations can't be negative...)
  real<lower=0, upper=1> theta;
  real<lower=0> n;
  real<lower=0, upper=1> m;
}
transformed parameters {
  vector[2] beta_pars;  // vector of parameters
  beta_pars[1] = n*m;
  beta_pars[2] = n*(1-m);
}
model {
  // This is where we write out the probability model, in very similar form to how we would using paper and pen

  // Define the priors
  n ~ gamma(s, r);
  m ~ beta(mu, nu);
  theta ~ beta(beta_pars[1], beta_pars[2]); // a= n*m, b=n*(1-m)
  
  // The likelihood
  X ~ binomial(N, theta);
}
generated quantities {
  // We will also generate posterior predictive draws (yhat) for each data point. These will be elaborated on below. 
  
  int y_sim;
  y_sim = binomial_rng(N, theta);
}"
```

```{r}
data_list <- list(N = nrow(chocolate), X = sum(chocolate$Good=="Yes"), nu = 2, mu = 4, s=4, r=0.2)
chocolate.hierarchical.fit <- stan(model_code = chocolate.hierarchical.model, data = data_list, warmup=2000, chains = 1, iter = 12000)
```

* Extract the parameters; Plot the posterior of $\theta$ and compare with the previous plots. 

```{r}
hierarchical.thetas <- extract(chocolate.hierarchical.fit, "theta", permuted=F)
theta.plot.data <- data.frame("theta"=c(hierarchical.thetas, thetas), "model" = rep(c("hierarchical", "simple"), each=10000))
ggplot(theta.plot.data, aes(x=theta)) + facet_grid(~model) + xlim(c(0.9, 1)) + geom_histogram(binwidth=0.002)
```

Note that the posterior from the hierarchical model is slightly wider. It's standard deviation is `r sd(hierarchical.thetas)` compared with `r sd(thetas)` for the simple model. This is reasonable, since we've added more uncertainity by adding next layer of randomness. 

* Plot the posterior distribution of parameters $n$ and $m$. Compare with the prior distributions.

```{r, include=F}
hier.n <- extract(chocolate.hierarchical.fit,"n", permuted=F)
hier.m <- extract(chocolate.hierarchical.fit,"m", permuted=F)
qplot(c(hier.n), xlim=c(0, 100), binwidth=1)
qplot(c(hier.m), binwidth=0.005, xlim=c(0, 1))
```

Note that the posterior distribution for $m$ is much different from the prior one. Our model has learned that chocolate is actually rather good. On the other hand, the posterior distribution for $n$ did not change that much. This is because in our data, there is not much information about how strongly $m$ should influence $\theta$. We'll get to that in the next section.


# Advanced Chocolate Modeling

So far, we've only inspected the overall proportion of tasty chocolate bars. The purpose of this task was to illustrate how to quantify the uncertainity in data. This is extremely important in medical applications, because we need to be certain that a drug has no adverse effects. But in the case of chocolate, it's rather useless. 

What's not useless, is using more advanced hierarchical models to compare the quality of chocolate of different companies. Such analyses are useful e.g. in market analyses. 

Let's assume that company $i$ has quality $\theta_i$. We define quality as the proportion of tasty chocolate bars produced by this company. Each company has it's own quality, but it's reasonable to assume that these qualities come from a single distribution. This distribution represents overall score of chocolates (or our prior beliefs about chocolate's tastiness regardless of the company it comes from). 

Such approach allows us to "exchange" information between companies. Suppose that we test two bars from a single company, and it turns out that both are tasty. Does this mean that the company's quality is 100\%? No! Just two bars is way too few to conclude that. However, we can conclude that this company's quality *might* be slightly higher than the average quality of the other companies. This would be reflected in the *posterior* distribution of this company's quality, $\theta_i$ - it would be slightly higher than the average posterior distribution of quality.

Let $X_i$ be the number of tasty chocolate bars produced by company $i$. We consider the following hierarchical model:

\begin{align}
X_i & \sim Bin(\theta_i, N_i) \\
\theta_i & \sim Beta(nm, n(1-m)) \\
n & \sim Gamma(s, r), \\
m & \sim Beta(\mu, \nu).
\end{align}

The only thing that's changed is that now we have multiple parameters $\theta_i$ (as many as the number of companies). Our data now consists of two vectors: first, encoding the numbers of tested chocolates of $i$-th company, and second, encoding the numbers of tasty bars. Now, the parameter $m$ corresponds to overall scoring of the chocolate bars, so it plays the role of parameter $\theta$ from the previous analyses. 

The parameter $n$ measures the influence of the prior scoring on the individual qualities. If $n$ is large (say, over 20), then the qualities $\theta_i$ are all similar to the overall quality $m$. This is because parameters $nm$ and $n(1-m)$ are large, and distributions of $\theta_i$ are very narrow (recall that $m$ is the mean value of $\theta_i$). On the other hand, if $n$ is small (say, around 5), then the companies differ greatly in their quality. Then, individual qualities $\theta_i$ can be very different than the overall quality $m$. This is because $nm$ and $n(1-m)$ are small, and the distributions of $\theta_i$ are very wide. Therefore, individual companies add a lot of variability to the overall distribution of tastiness given by $m$. 

We will assume a fairly weak prior (i.e. one that does not give much information). For this, we will let $n$ be equal to $10 \pm 10$, which corresponds to $s = 1$ and $r = 0.1$. This means that our prior information is worth on average $10 \pm 10$ chocolate bars, and that it is most probably not worth more than $50$ bars (because the density of Gamma distribution is small for values larger than 50). Next, we will assume that the mean quality $m$ is most probably between $0.6$ and $1$, and most probably not lower than $0.5$. These conditions are satisfied by a prior on $m$ with $\mu = 10$ and $\nu=2$. Feel free to adjust those parameters so that they reflect your own beliefs about chocolate ratings. 

**Your tasks:**

* Create a data frame consisting of three columns: company name, tested bars, tasty bars.

    * Remember to use `tapply(DATA, INDEX, FUN)` to apply function `FUN` to subvectors of `DATA` determined by values in vector `INDEX` 
    * Tapply returns a named vector, in which names correspond to values of `INDEX` variable. Those names can be retrieved by function `names(x)`. You can use this feature to obtain the vector of company names.
    
```{r, include=F}
choc.nbs <- tapply(chocolate$Good, chocolate[,1], length)
choc.ok <- tapply(chocolate$Good, chocolate[,1], function(x) sum(x=='Yes'))
company.data <- data.frame("Brand" = names(choc.nbs), "Number"=choc.nbs, "Tasty"=choc.ok)
```

```{r}
head(company.data)
```

* Implement the hierarchical model in STAN. 

    * The model should accept number of companies, a vector of numbers of bars tested per brand, numbers of tasty bars, and hyperparameters $s, r, \nu$ and $\mu$.
    * DO NOT generate the predicted distribution of data in this example (the variable `y_sim`).
    * The default type of vectors in STAN is real. To get an integer vector of length `L`, write `int X[L]`.
    * Syntax of for loop in STAN is the same as in R, but without braces (e.g. `for(i in 1:L) X[i] ~ binomial(N[i], theta[i]);`).
    
* Run the analysis. Extract the parameters. Use 12 000 iterations, including 10 000 warmup iterations.
    
```{r, include=F}
company.model <- "
data {
  // In this section, we define the data that must be passed to Stan (from whichever environment you are using)

  int L; // number of companies
  int N[L]; // numbers of bars tested
  int X[L]; // number of tasty bars
  real s;
  real r; // hyperparameter of n 
  real mu;
  real nu; // hyperparameters of m
}
parameters {
  // Define the parameters that we will estimate, as well as any restrictions on the parameter values (e.g. standard deviations can't be negative...)
  vector<lower=0, upper=1>[L] theta;
  real<lower=0> n;
  real<lower=0, upper=1> m;
}
transformed parameters {
  vector[2] beta_pars;  // vector of parameters
  beta_pars[1] = n*m;
  beta_pars[2] = n*(1-m);
}
model {
  // This is where we write out the probability model, in very similar form to how we would using paper and pen

  // Define the priors
  n ~ gamma(s, r);
  m ~ beta(mu, nu);
  for(i in 1:L)
    theta[i] ~ beta(beta_pars[1], beta_pars[2]); // a= n*m, b=n*(1-m)
  
  // The likelihood
  for(i in 1:L)
    X[i] ~ binomial(N[i], theta[i]);
}
"
```


```{r}
data_list <- list(L=nrow(company.data), N = company.data$Number, X = company.data$Tasty, nu = 2, mu = 10, s=1, r=0.1)
company.fit <- stan(model_code = company.model, data = data_list, warmup=10000, chains = 1, iter = 15000, control = list(max_treedepth = 25))
```

```{r}
posterior.n <- extract(company.fit, "n", permuted=F)
posterior.m <- extract(company.fit, "m", permuted=F)
posterior.thetas <- extract(company.fit, "theta", permuted=F)
```


* Plot the posterior distributions of $n$ and $m$, and the dependence between those two parameters. How do they differ from the prior ones?

```{r}
qplot(c(posterior.n), c(posterior.m))
```

* Plot the distributions of several chosen $\theta_i$ using `facet_wrap`. Interpret the results. 

```{r}
my.companies <- data.frame(posterior.thetas[,,c(8, 120, 232, 398)])
my.companies <- melt(my.companies)
ggplot(data=my.companies, aes(x=value)) + facet_wrap(~variable) + geom_histogram(binwidth=0.002)
```

Let's choose company number 232, "Machu Picchu Trading Co.", to illustrate an important phenomenon. You can see that the plot of this company's quality looks different - it clearly seems that this one is not doing as well as the others. If you inspect the data, you'll see that two bars of this manufacturer were tested, and only one was tasty. However, note that the value 0.5 is not very probable in the posterior distribution. The so-called "90% Highest Density Interval" -- the shortest interval encompassing 90% of the probability -- is equal to $[0.82, 0.99]$. This means that in 90% of iterations, the quality belonged to this interval. Of course, it is possible that the quality is below 0.5 -- it's very unlikely, but possible. 

This illustrates the phenomenon of information transfer between groups. We know that most of the bars are tasty. Therefore, even though only half of the *tested* bars of Machu Piccu Trading Co. are tasty, it's unlikely that only half of *all* the bars of this company are tasty. After all, we've only tested two of them, and there may be hundreds of other tasty brands from this company. All we know is that this company may have lower quality than the other companies. How much lower? This is expressed by the 90% HDI - we can say that we're 90% certain that the quality is between 0.82 and 1. 

Compare this with the company number 8 ("Altus aka Cao Artisan"). In this case, based on the 90% HDI, we're 90% certain that the quality is between 0.95 and 1. For this company, we've tested 10 bars and all were tasty. The same princile of information transfer applies: Even though all the bars were tasty, we're not 100% sure that the quality is equal to 1, and all the brands of this company are equally good. In fact, it's very unlikely that this company has never ever produced an untasty bar. But we're 90% sure that the quality is very high -- above 0.95.

We won't cover the HDI in more detail now, but here's a function to compute it if you're interested:

```{r}
HDIofMCMC = function( sampleVec , credMass=0.9 ) {
  # Computes highest density interval from a sample of representative values,
  #   estimated as shortest credible interval.
  # Arguments:
  #   sampleVec
  #     is a vector of representative values from a probability distribution.
  #   credMass
  #     is a scalar between 0 and 1, indicating the mass within the credible
  #     interval that is to be estimated.
  # Value:
  #   HDIlim is a vector containing the limits of the HDI
  sortedPts = sort( sampleVec )
  ciIdxInc = ceiling( credMass * length( sortedPts ) )
  nCIs = length( sortedPts ) - ciIdxInc
  ciWidth = rep( 0 , nCIs )
  for ( i in 1:nCIs ) {
    ciWidth[ i ] = sortedPts[ i + ciIdxInc ] - sortedPts[ i ]
  }
  HDImin = sortedPts[ which.min( ciWidth ) ]
  HDImax = sortedPts[ which.min( ciWidth ) + ciIdxInc ]
  HDIlim = c( HDImin , HDImax )
  return( HDIlim )
}
```




