---
title: "SISG Bayes: Session 5 FTO/Stan example"
author: "Ken Rice"
date: "`r Sys.Date()`"
output: pdf_document
---

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

*Before you start:* please check the course website for instructions: Windows users will need the RTools bundle, and the rstan package download is non-standard. 

Once those steps are done, load the data:

```{r}
# install.packages("rstan")
setwd("C:/Users/kenrice/OneDrive - UW/Desktop/SISGBayes/exercises")
load("yX_FTO.Rdata")
y <- yX$y
X <- yX$X
n <- nrow(X)
p <- ncol(X)
```
Write a file containing a description of the model in Stan's language, it has elements of C, R and BUGS/JAGS:

```{r}
cat(file="FTOexample.stan", "
data {
  int n; //the number of observations
  int p; //the number of columns in the model matrix
  real y[n]; //the response
  matrix[n,p] X; //the model matrix
  real g; // Zellner scale factor
  vector[p] mu; // Zellner prior mean (all zeros)
  matrix[p,p] XtXinv; // information matrix
}
parameters {
  vector[p] beta; //the regression parameters
  real invsigma2; //the precision, a.k.a. inverse-variance
}
transformed parameters {
  vector[n] linpred;
  cov_matrix[p] Sigma;
  real sigma;
  linpred = X*beta;
  sigma = 1/sqrt(invsigma2);
  for (j in 1:p){
    for (k in 1:p){
      Sigma[j,k] = g*sigma^2*XtXinv[j,k];
      } } }
model {  
  beta ~ multi_normal(mu, Sigma);
  y ~ normal(linpred, sigma);
  invsigma2 ~ gamma(0.5, 1.839); // we took nu0=1, sigma0=1.91
}
")
```

From within R, we set Stan going -- first compiling the model code, to write code that does the MCMC sampling, then running that code and collecting its output in an R object:

```{r}
library("rstan")
stan2 <- stan(file = "FTOexample.stan", 
data = list(n=n,p=p, y=y, X=X, g=n, mu=rep(0,p), XtXinv=solve(crossprod(X)) ),
iter = 100000, chains = 1, pars=c("beta","sigma"))
```
Some summaries of the posterior samples -- using packaged code instead of `by hand'.

```{r}
print(stan2)
plot(stan2)
traceplot(stan2)
```