Observed data (simulated, of course!)

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.

sig2.mle = mean(x^2)
sig2.mle
## [1] 0.2598108

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\):

#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

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)
}
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)\)

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
## [1] 0.2518762
sig2.mode = b/(a+1)
sig2.mode
## [1] 0.2598108

Riemann’s integration for \(Pr(\sigma^2>0.25|x_1,\ldots,x_n)\)

(pinvgamma(u,a,b)-pinvgamma(0.25,a,b))/pinvgamma(u,a,b)
## [1] 0.6760581

Posterior quantiles (via inverse CDF transformation)

alpha = c(0.975,0.5,0.025)
ci = qinvgamma((1-alpha)*pinvgamma(u,a,b),a,b)
ci
## [1] 0.1411906 0.2986232 0.4843162
sig2.median = ci[2]

# Estimates: MLE & posterior mean, mode and median
c(sig2.mle,sig2.mean,sig2.mode,sig2.median)
## [1] 0.2598108 0.2518762 0.2598108 0.2986232

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\).