---
title: "HW05 Solutions"
author: "BJ"
date: '2022-03-08'
output: pdf_document
---

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

Install and verify the software package Stan and the library rstan, on your personal computer. Instructions for doing this are appended at the end of this hw assignment.  Then run the following short example to verify that rstan is intalled and working correctly (note that the text string assigned to the variable ``model'' runs over several lines; that's perfectly OK):
```{r warning=F, message=F}

set.seed(1234567890) ## to make a reproducible simulation result 

library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
library(ggplot2); theme_set(theme_bw())
library(bayesplot)
```


```{r}
text_model <- "
   data {
     real y_mean;
   } 
   parameters {
     real y;
   } 
   model {
     y ~ normal(y_mean,1);
   }"

```
This Stan program is not doing anything Bayesian, it's just making draws from a normal distribution with mean=0 and sd=1 (fit) or mean 5 and sd=1 (fit2).

## Problem 1(a).

Turn in the output for both \verb|print(fit)| (or just \verb|fit|), and \verb|summary(fit)|.  Write a sentence or two explaining the difference between the two outputs (besides a different number of decimal places printed!).



```{r}
compiled_model <- stan(model_code=text_model,data=list(y_mean=0),chains=0)
fit <- stan(fit=compiled_model,data=list(y_mean=0))
fit2 <- stan(fit=compiled_model,data=list(y_mean=5))
```


```{r}
## fit or print(fit)
fit
## summary(fit)
summary(fit)
```

\it
Here's a comparison of the two outputs:
\begin{itemize}

\item \verb|print(fit)| just gives a summary of the estimated parameters, \verb|se_mean| from the Monte Carlo sumualtion, the posterior \verb|sd|, quantiles or the posterior distribution, and \verb|n_eff| \& \verb|Rhat| values, for the "kept" part of all the chains, merged together.

\item \verb|summary(fit)| gives the same overall summary, as a two-dimensional array, in the \verb|$summary| component, and then in \verb|$c_summary| it provides a 3-dimensional array containing summaries for each individual Markov chain that was run (but without the \verb|n_eff| \& \verb|Rhat| values). 

\item[] (The two-dimensional array in \verb|summary(fit)$summary| is convenient for extracting mean, sd, quantiles, etc. to make CI's. etc.)

\end{itemize}

\rm

## Problem 1(b).

Use commands like those in the appendix (from \verb|library(bayesplot)|) of lecture 13 to produce and turn in, for the fitted model \verb|fit|:
\begin{itemize}
\item Trace plots of \verb|y| (the variable being simulated) and \verb|lp__| (the value of the log-posterior (in this case, not a posterior but just the normal density) at each simulated value of \verb|y|). 
\item Autocorrelation plots of \verb|y| and \verb|lp__|.
\item A histogram of just \verb|y|. 
\end{itemize}

```{r fig.height=5}
library(gridExtra)
g1 <- mcmc_trace(fit)
g2 <- mcmc_acf(fit)
g3 <- mcmc_areas(fit,pars="y")
grid.arrange(arrangeGrob(g1,g2,ncol=2),g3,heights=c(2.5/4,1.4/4),ncol=1)
```


## Problem 1(c).

Use the \verb|extract| function to extract simulated values for \verb|y| from \verb|fit| and make a qqplot of the values against the Normal distribution.  Is it plausible that the simulated values for \verb|y| are normally distributed?

```{r fig.height=2.75}
sim <- extract(fit)

ggplot(as.data.frame(sim),aes(sample=y)) +
  geom_qq() +
  geom_qq_line() +
  xlab("Normal Quantiles") + 
  ylab("y (from \"fit\")")
```
\it

Yes, given this plot, it's pretty plausibly that the y's are normally distributed. There's just a hint of a deviation from the normal in the tails (where there's very little sample size, so it's hard to say these are "significant" deviations), but overall the data follow the normal distribution well.

\rm

## Problem 1(d).

Print out and turn in output from \verb|print(fit)| and \verb|print(fit2)|.  Write a sentence noting any similarities or differences between them.


```{r}
fit
fit2
```

\it

Everything is pretty similar for the two outputs, except of course the estimated mean (as either a posterior mean or a posterior median) is $\approx$ 0 for \verb|fit| and it is $\approx$ 5 for \verb|fit2|.  The posterior quantiles are similarly shifted by about 5 units.

\rm

## Problem 1(e).

Use the \verb|extract| function to extract simulated values for \verb|y| from \verb|fit2| and make a qqplot of the values against the Normal distribution.  Is it plausible that the simulated values for \verb|y| are normally distributed?
```{r fig.height=2.75}
sim <- extract(fit2)

ggplot(as.data.frame(sim),aes(sample=y)) +
  geom_qq() +
  geom_qq_line() +
  xlab("Normal Quantiles") + 
  ylab("y (from \"fit2\")")
```

\it

Again this qq plot is pretty convincing of normality.  There's just a hint of a deviation from the normal in the tails (where there's very little sample size, so it's hard to say these are "significant" deviations).  Note that the y-intercept for the fitted line is now about 5, corresponding to the fact that we are simulating from a $N(5,1)$ instead of a $N(0,1)$ this time.

\rm

# Problem 2.

Problem 2(c)(iv) on hw \#04 asked you to generate a point estimate and 95\% CI for the model
\begin{eqnarray*}
    y_i &\sim & Poiss(\lambda), ~i=1,2, \ldots, n \\
    \lambda & \sim & Gamma(\alpha,\beta)
\end{eqnarray*}
with $\alpha=1$ and $\beta=1$, for the $n=16$ data values
\[
	3,  7,  1,  7,  6,  1,  5,  4,  6,  2, 11,  3,  1,  6,  4,  4
\]

## Problem 2(a).

Write a stan program that specifies the likelihood and prior for this model, and use it to simulate draws from the posterior.  You can use the default number of chains (4) and samples in each chain (2000).  Include in your solutions:
\begin{itemize}
\item Your stan program
\item your R code
\item the output from \verb|print()|'ing fitted stan simulation.
\end{itemize}
\textit{Hint:}  You will find the documentation at 

\url{https://mc-stan.org/docs/2_29/functions-reference/index.html}

useful, especially the sections on Discrete and Continuous Distributions.  You may also find this documentation useful

\url{https://mc-stan.org/docs/2_29/stan-users-guide/index.html}

\url{https://mc-stan.org/docs/2_29/reference-manual/index.html}

and, of course, you should imitate examples given in class to save yourself time, wherever possible.


```{r}
## Either of the following Stan programs works. The only difference is whether you 
## let stan "vectorize" the part of the code that establishes the poisson distribution 
## for each observation y[i], or you do it yourself with a "for" loop...

text_model <- "
data {
  int n;
  int y[n];
}
parameters {
  real lambda;
}
model {
  y ~ poisson(lambda);   // likelihood
  lambda ~ gamma(1,1);   // prior
}"

text_model <- "
data {
  int n;
  int y[n];
}
parameters {
  real lambda;
}
model {
  for (i in 1:n) {       // likelihood
    y[i] ~ poisson(lambda);   
  }
  lambda ~ gamma(1,1);   // prior
}"

## Then we compile and run the model as follows:

compiled_model <- stan(model_code=text_model,
                       data=list(n=2,y=c(1,2)),
                       chains=0)

fit <- stan(fit=compiled_model,
            data=list(y=c(3,7,1,7,6,1,5,4,6,2,11,3,1,6,4,4),
                      n=16))

print(fit)
```


## Problem 2(b).

Use functions in \verb|library(bayesplot)| to produce and turn in graphs of
\begin{itemize}
\item trace plots
\item acf plots
\item rhat plots
\item neff (ratio) plots
\end{itemize}
for your fitted simulation.

```{r}
g1 <- mcmc_trace(fit)
g2 <- mcmc_acf(fit)
g3 <- mcmc_rhat(rhat(fit)) + yaxis_text()
g4 <- mcmc_neff(neff_ratio(fit)) + yaxis_text()

grid.arrange(g1,g2,g3,g4,ncol=2)
```


## Problem 2(c).

Use functions in \verb|library(bayesplot)| to produce and turn in graphs of 
\begin{itemize}
\item A 95\% CI for $\lambda$
\item A density plot for $\lambda$, with the 95\% CI shaded in underneath it.
\end{itemize}
Compare your answer here with the answers you got on hw04 problem \#2.

```{r}
g1 <- mcmc_intervals(fit,pars="lambda",prob_outer=0.95)
g2 <- mcmc_areas(fit,pars="lambda",prob_outer=0.95)
grid.arrange(g1,g2,ncol=1)
```
\it

Here are the intervals and point estimates from HW04:
\begin{verbatim}
## sample avg and sd:
round(c(xbar + c(-2,0,2)*stdev/sqrt(n)),2)
## [1] 3.08 4.44 5.79
## MLE:
round(c(lambda_hat + c(-2,0,2)*se_lambda_hat),2)
## [1] 3.38 4.44 5.49
## glm:
round(exp(c(beta_hat + c(-2,0,2)*se_beta_hat)),2)
## [1] 3.50 4.44 5.63
## gamma posterior, equal-tailed:
round(qgamma(c(0.025,0.5,0.975),astar,bstar),2)
## [1] 3.31 4.22 5.27
## gamma posterior, wald interval:
round(c(gamma_mean + c(-2,0,2)*gamma_sd),2)
## [1] 3.24 4.24 5.23
\end{verbatim}
The interval from our simulation,
```{r}
round(summary(fit)$summary[1,c(4,6,8)],2)
```
is very similar to the interval we calculated using quantiles of the Gamma distribution.  (This isn't at all suprising, since the posterior distribution we are simulating from with stan is exactly the same Gamma distribution we calculated for hw04).

\rm


## Problem 3.

In the files area of Canvas there is a folder ``presentation project ideas''. 

## Problem 3(a).

Read the ``overview.pdf'' file in that folder, and browse through the subfolders there\footnote{If there is only a readme.txt file there now, it means I am still organizing and uploading content. It should all be there by Wednesday evening.}.  If you might have a project of your own to present, think about that too.  Turn in your first, second and third choices for projects (these are tentative, not final, choices!).  Also, indicate tentative team members if you are organizing a team.

\it

If you list 3 ideas for presentation projects, and list team members (if you are part of a team), you will get a full 10 pts for this part.

\rm

## Problem 3(b).

Sign up for a time to meet me in the week after spring break to talk about your presentation project. (You don't have to turn anything in for this, but it will be counted as part of your project grade.)

\it

You get 10 automatic free points on the hw for this part; you don't have to turn anything in for these 10 pts!

I will separately give you a score in the "presentation projects" part of the gradebook for our meeting to discuss which presentation project you will actually do.

\rm

\clearpage 

\begin{center}
\Large \textbf{ Installing Stan  for R}
\end{center}

\begin{enumerate}

\item Go to \url{https://mc-stan.org/users/interfaces/rstan} \\ (can also find by
   googling ``rstan'').

\item Click on the ``RStan Quick Start Guide'' under ``Download and Get
   Started''

   (\textit{Note:} We will be installing Stan to run in R.  There are also
   installations for python, matlab, julia, stata, mathematica, scala,
   etc\ldots{}  If you are interested, see \url{https://mc-stan.org/} for
   details.)

\item Follow the instructions for configuring the \verb|C++| Toolchain on your
   computer
   
   \begin{itemize}

   \item This is the trickiest part of installing Stan in R on your
     computer.

   \item There are separate instructions for Windows, Mac, \& Linux

   \item Make sure you are using R version 4.0 or higher.  You may want to
     update Rstudio to the current version as well (you need at least
     version 1.2.5042).

   \item \textbf{ MAKE SURE YOU ARE} following the instructions for your version of
     R (version 4.0 or higher) not an earlier version of R such as
     version 3.6, etc.
     
   \end{itemize}
     
\item Follow the instructions under ``Installation of RStan'' and ``Verify
   Installation''
   
   \begin{itemize}

   \item The most important parts are the ``install.packages'' command, to
     actually do the installation, and the ``example'' command to verify
     that your installation is good.  If you can get these two
     commands to run successfully, you are golden.
     
     \begin{itemize}

     \item Note: the ``example'' command will take from several seconds to
        up to a minute or so to run.  It will also generate a ton of
        output that may not make much sense right now\ldots{} Be patient\ldots{}

     \end{itemize}

   \item Please note the other instructions under these two headings,
      which handle some common special cases that you may need to deal
      with on your laptop.

   \end{itemize}

\item After you have verified that your installation is good, feel free
   to browse through the longer ``How to Use RStan'' section.  As you
   get more accustomed to Stan, you'll see that the Stan modeling
   language is a kind of extension of R's usual programming language
   -- much more general and flexibe than \verb|lmer()|'s modeling language.

\item Check out \url{https://mc-stan.org/} for background info on stan.

\end{enumerate}

