Homework Set 16, Solutions      11/21/2016
------------------------------------------

> Ydat16 = scan("http://www.math.umd.edu/~slud/s705/Homework/HW16dat.txt")

## posterior is mu ~ N( 16n Ybar/(16n+1), 16/(16n+1) )

set.seed(355)
mu.st = rnorm(1000, (480/481)*mean(Ydat16), sqrt(16/481))

> hist(mu.st,nclass=32)    ### doesn't actually look very normally distributed

Yst.arr = array(rnorm(30*1000, rep(mu.st,each=30), 1), c(30,1000))   
    ### 30 x 1000, with predictive data-batches in columns

> summary(apply(Yst.arr,2,var))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.3583  0.8081  0.9533  0.9780  1.1480  1.9560 

### The variances of the batches look quite like 1

> summary(apply(Yst.arr,2,mean))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  6.344   7.034   7.213   7.204   7.379   7.989

### and the means are all roughly centered around Ybar = 7.21757

### So it will not be a surprise that the means of the predictive Y.star batches
###   surround the mean of the original dataset.

## How about 2nd moments ?

ymom2 = apply(Yst.arr^2,2,mean)

> c(y2ndmom = mean(Ydat16^2), CI.star = mean(ymom2) + 1.96*c(-1,1)*sd(ymom2)/sqrt(1000))
 y2ndmom CI.star1 CI.star2 
70.19272 52.67969 53.13156 

### We could also do it with 3rd moments, which will be even more extreme

ymom3 = apply(Yst.arr^3,2,mean)

> c(y3rdmom = mean(Ydat16^3), CI.star = mean(ymom3) + 1.96*c(-1,1)*sd(ymom3)/sqrt(1000))
 y3rdmom CI.star1 CI.star2 
853.6438 393.1836 398.1879 

#### Let's finish off with posterior predictive checking of proportion  > 8

ygt8 = apply(Yst.arr>8, 2, mean)

>  c(ygt8 = mean(Ydat16>8), CI.star = mean(ygt8) + 1.96*c(-1,1)*sd(ygt8)/sqrt(1000))
     ygt8  CI.star1  CI.star2 
0.3666667 0.2076438 0.2188895

### So we can see pretty convincingly in this case that posterior predictive checking
###   says the fit of this model is no good !!