---
title: "Gibbs Sampling and Getting What You Paid For"
author: "Carolyn J. Anderson"
date: "9/2/2021"
output: 
   html_document:
     toc: TRUE
---
# Introduction
Goal for edsp690BAY: Find the mean and variance of the state average total SAT total score (i.e., verbal + math) using Gibbs sampling.

This gives you a start and we'll go over this and finish in class.

In this practice problem, we are using data from  "Getting What Paid For: Debate Over Equity in Public School Expenditures" which contains state data measures on cost of education and SAT scores.  These data are from 1994-1995.  We will use this data later when talk about linear regression, but for now, we'll just estimate the mean and assume that we know the variances.

The data are from Deborah Guber and the full reference is Guber, D. (1999). "Getting What You Pay For: The Debate Over Equity in Public School Expenditures". Journal of Statistics Education.  For more information about the data see http://jse.amstat.org/datasets/sat.txt.


# Set-up



## Tools
We'll use the following libraries
```{r setup, message=FALSE, warning=FALSE, include=TRUE}
library(coda)
library(rjags)
library(runjags)
library(HDInterval)  # for high density interval
```

You might also want to use functions from Kruschke, so I put these in as source code.
```{r}
source("C:/Users/cja/Dropbox/edps 590BAY/Lectures/5 MCMC Gibbs/DBDA2E-utilities.txt", encoding = 'UTF-8')
```

## Data

You can find the data on the course web-site.  You need to change the path to where ever you put the data.

```{r}
df <- read.table("C:/Users/cja/Dropbox/edps 590BAY/Lectures/5 MCMC Gibbs/getting_what_paid_for_data.txt", header=TRUE)
names(df)
head(df)
```

# Plot the Data

..and overlay the normal distribution.  There are a number of ways to do this and what's below keeps things in terms of frequencies.  The first line sets up the histogram and saves it as "h".  The next two lines get data for horizontal and vertical axis.  The third command puts density values of vertical fit values into scale of frequencies and uses information from the historgra.  The "lines" command puts in the normal curve.  

```{r}
h =hist(df$ave_tot, breaks = 10, density = 10,
            col = "darkblue", 
            xlab = "Average Total SAT", 
            ylab = "Frequency",
            main="Normal curve Over histogram")
xfit <- seq(min(df$ave_tot), max(df$ave_tot), length = 50) 
yfit <- dnorm(xfit, mean = mean(df$ave_to), sd = sd(df$ave_tot)) 
yfit <- yfit * diff(h$mids[1:2]) * length(df$ave_tot) 
lines(xfit, yfit, col = "black", lwd = 2)
```

# Sample Statistics

We need the sample size and observed sample mean.  We will pretend that we know what the variance equals and use the sample variance for this.

```{r}
n <- nrow(df)
ybar <- mean(df$ave_tot)
std <- sd(df$ave_tot)
s2 <- var(df$ave_tot)


descriptive <- c(n, ybar, std, s2)
names(descriptive) <- c("N","Ybar","std","s^2")
descriptive
```

# Using rjags

## model

## inital/starting values

## Run rjags

## Assess convergence

# Use runjags
