---
title: "Estimating Mean and Variance for Normal"
author: "Carolyn J. Anderson"
date: "9/3/2021"
output: 
   html_document:
       toc: yes
---
# 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 and variance 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

This is pretty much the same as we did for when estimating mean but the variance was known.

## 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)
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
```

# Flat (diffuse) and non-informative Priors

We'll start with flat priors reflecting no knowledge and ignorance about SAT scores and these data. 

```{r}
mu.o <- 0          #  mean<<possible minimum of 400 
tau2.o <- 1000000  # Very flat, N(0, 1000000)
kappa.o <- 1
nu.o <- 1
sigma2.o <- 1000000
```


And the posterior distribution of the mean is normal with mean mu.n and variance tau.n; that is
$$ \bar{y} \sim N(\mu_n, \tau_n^2) \qquad \mbox{and} \qquad 1/\sigma^2 \sim \mbox{Gamma}(\nu_{n}/2, \  (\nu_{n}/2)\sigma^2_{n})$$
will be our posterior distribution.

```{r}
# Posterior sample size
(kappa.n <- kappa.o + n)

# Posterior nu
(nu.n <- nu.o + n)

# Posterior mean (mu.n)
(mu.n <- (kappa.o*mu.o + n*ybar)/(kappa.o + n) )
  
(sigma2.n <- (nu.o*sigma2.o + (n-1)*s2 + (kappa.o*n/kappa.n)*(ybar-mu.o)**2)/nu.n)

(s2_theta.n <- sigma2.n/nu.n)

(tau2.n <- 1/(kappa.o/tau2.o + (nu.o + n)/sigma2.n))
```


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

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

The results don't seem right and the problem is the poor priors.


# Priors Using Some Information About SAT scores


I know that SAT total scores range from 400 to 1600 (i.e., 200 min and 800 max on each of the two tests). 
A score of 500 on quantitative and verbal tests might be reasonable averages for each test, so I set mu.o to 1000. I won't change any of the other priors
```{r}
mu.o <-1000   
tau2.o <- 1000000  # Very flat, N(1000, 1000000) 
kappa.o <- 1
nu.o <- 1
sigma2.o <- 10000
```

```{r}
# Posterior sample size
kappa.n <- kappa.o + n

# Posterior nu
nu.n <- nu.o + n

# Posterior mean (mu.n)
mu.n <- (kappa.o*mu.o + n*ybar)/(kappa.o + n) 
  
sigma2.n <- (nu.o*sigma2.o + (n-1)*s2 + (kappa.o*n/kappa.n)*(ybar-mu.o)**2)/nu.n
```

Two ways of computing the variance of the mean, lead to about same result
```{r}
s2_theta.n <- sigma2.n/nu.n

tau2.n <- 1/(kappa.o/tau2.o + (nu.o + n)/sigma2.n)
```

These results look much more reasonable and comparable to sample statistics.  As a check, I used GIBBS (rjags) and get similar results.

Putting this all together,
$$ \mu_n \sim N(`r round(mu.n,2)`,`r round(tau2.n,2)`)\quad \sigma^2_n \sim \mbox{IGamma}(`r nu.n`, `r round(sigma2.n, 2)`)\quad \mbox{ or }\quad   1/\sigma_n^2 \sim \mbox{Gamma}(`r nu.n/2`, `r nu.n*sigma2.n/2`)$$
Note standard deviation of $\mu_n$ is $'r \sqrt{tau2.n}'$

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

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


# Monte Carlo 

We first draw random samples for precision (i.e., 1/$\sigma^2$). Then we draw random  values from the normal using the simulated values of $\sigma^2$.   Taking 1000 random draws from the posterior

```{r}
sigma2.mc <- 1/ rgamma(1000, nu.n/2, nu.n*sigma2.n/2)
tau2.mc <- sigma2.mc/(kappa.o + n)
mu.mc <- rnorm(1000, mu.n, sqrt(tau2.mc))
```

Compare the following with the analytic estimates of the mean, variance and precision,
```{r}
summary(mu.mc)
summary(tau2.mc)
summary(sigma2.mc)

sd(mu.mc)
var(mu.mc)
```

Pretty close. 
