---
title: "Estimating Mean from Normal with Variance Known"
author: "Carolyn J. Anderson"
date: "9/2/2021"
output: 
   html_document:
     toc: TRUE
---
# Introduction

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.

Goal for edsp690BAY: Find the mean of the state average total SAT total score (i.e., verbal + math).  This gives you a start and we'll go over this and finish in class.


# Set-up

## Tools
Library that I used
```{r setup, message=FALSE, warning=FALSE, include=TRUE}
library(HDInterval)  # for high density interval
```

## 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("D:/Dropbox/edps 590BAY/Lectures/3 Normal/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))
(sigma2 <- var(df$ave_tot))
```


# Flat Priors

We'll start with flat proiors refelcting no knowledge about SAT scores. 

```{r}
mu.o <- 0     
tau.o2 <- 1000000  # Very flat, N(0, 1000000):
```

And the posterior distribution of the mean is normal with mean mu.n and varianace tau.n; that is
$$ \bar{y} \sim N(\mu_n, \tau_n^2)$$
will be our posterior distribution.

```{r}
# Posterior Mean 
(mu.n <- ((1/tau.o2)*mu.o + (n/sigma2)*ybar)/(1/tau.o2 + n/sigma2))

# Posertior variance
(tau.n2 <-  1/ ( 1/tau.o2 + n/sigma2) )
```

Intervals estimates
```{r}
# Credible interval
(pq <- qnorm(c(.025, .975), mu.n, sqrt(tau.n2)))

# High Density interval
x <- rnorm(100000, mu.n, sqrt(tau.n2))
library(HDInterval)
hdi(x, credMass=.95)
```


# Use What you Know about SAT Scores 

For this,


* Use what you know about SAT scores and change the prior accordingly. 
* Compute credible and high density intervals.
* Compare to your results with non-informative priors (i.e., ignorance)

# Try out other Priors

To see what effect they have on results.


# For Shadow Future Analsyses

```{r}
names(df)
cor(df$exp_pp, df$ave_tot)

plot(df$salary, df$ave_tot, type="n", 
     main="Getting what paid for?  r(expendature/pupil, ave_total_SAT) = -.38",
     xlab="Expenditure per Student",
     ylab="Average Total SAT Score",
     xlim=c(20,55))
text(df$salary, df$ave_tot, df$state_abv)
```