---
title: "Posterior inference for the variance in an IID Normal problem"
author: "Hedibert Freitas Lopes"
date: "5/4/2018"
output:
# word_document: default
# pdf_document: default
html_document: default
---

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

## Observed data (simulated, of course!)

```{r}
set.seed(1234)
n    = 10
sig2 = 0.25
x    = rnorm(n,0,sqrt(sig2))
```

## Statistical model
We assume that, conditional on $\sigma^2$, the observations $x_1,\ldots,x_n$ are iid Gaussian with mean zero and variance $\sigma^2$.  We are interested in cases where $\sigma^2<u$, for a given upperbound $u>0$.  In this case, it is easy to see that the MLE of $\sigma^2$ is 
$$
{\hat \sigma}_{MLE}^2 = \frac{1}{n}\sum_{i=1}^n x_i^2,
$$
for ${\hat \sigma}_{MLE}^2<u$, or $u$ otherwise.

```{r}
sig2.mle = mean(x^2)
sig2.mle
```

## Bayesian estimation
We will use a uniform prior for $\sigma^2$, i.e., $U(0,u)$.  In what follows, we use $u=0.5$.
In this situation the posterior distributionof $\sigma^2$ is a truncated inverse gamma:
$$
\sigma^2 | x_1,\ldots,x_n \in IG\left(\frac{n-2}{2},\frac{\sum_{i=1}^n x_i^2}{2}\right)
1_{(0,u)}(\sigma^2).
$$
The posterior density is, therefore,
$$
p(\sigma^2 | x_1,\ldots,x_n) = \frac{f_{IG}(\sigma^2;a,b)}{F_{IG}(u;a,b)},
$$
where, $a=(n-2)/2$, $b=\sum_{i=1}^n x_i^2/2$, $f_{IG}(\cdot;a,b)$ and $F_{IG}(\cdot;a,b)$ are, respectively, the pdf and cdf of the inverse gamma with parameters $a$ and $b$.

# A few R functions to evaluate the truncated inverse gamma posterior of $\sigma^2$:
```{r}
#install.packages("invgamma")
library(invgamma)
dinvgamma.trunc=function(sig2,a,b,u){if(sig2<u){dinvgamma(sig2,a,b)/pinvgamma(u,a,b)}else{0}}
```

We are now ready to plot the posterior density of $\sigma^2$:
# Plot of the posterior density
```{r}
u = 0.5
a = (n-2)/2
b = sum(x^2)/2
N = 1000
sig2s = seq(0.0001,2*u,length=N)
den = rep(0,N)
den.trunc = rep(0,N)
for (i in 1:N){
 den[i] = dinvgamma(sig2s[i],a,b)
 den.trunc[i] = dinvgamma.trunc(sig2s[i],a,b,u)
}
```

```{r fig.width=8, fig.height=5}
par(mfrow=c(1,1))
plot(sig2s,den.trunc,type="l",lwd=3,xlab=expression(sigma^2),ylab="Posterior density")
lines(sig2s,den,type="l",col=2,lwd=3)
abline(v=sig2,col=4,lwd=3,lty=3)
legend("topright",legend=c("True sig2","IG","Truncated IG"),col=c(4,2,1),lty=c(3,1,1),lwd=3)
```

# Riemann's integration for $E(\sigma^2|x_1,\ldots,x_n)$
```{r}
num = 0.0
for (i in 1:N)
  num = num + sig2s[i]*dinvgamma(sig2s[i],a,b,u)*u/N
den = pinvgamma(u,a,b)
sig2.mean = num/den
sig2.mean

sig2.mode = b/(a+1)
sig2.mode

```

# Riemann's integration for $Pr(\sigma^2>0.25|x_1,\ldots,x_n)$
```{r}
(pinvgamma(u,a,b)-pinvgamma(0.25,a,b))/pinvgamma(u,a,b)
```

# Posterior quantiles (via inverse CDF transformation)
```{r}
alpha = c(0.975,0.5,0.025)
ci = qinvgamma((1-alpha)*pinvgamma(u,a,b),a,b)
ci
sig2.median = ci[2]

# Estimates: MLE & posterior mean, mode and median
c(sig2.mle,sig2.mean,sig2.mode,sig2.median)
```

# Your turn
You can now repeat the above exercise by changing the sample size, $n$, the prior upperbound,$u$, and the true value of $\sigma^2$.