rm(list=ls())

#-----------------------------------------------------#
#---   Classical and Bayesian inference (AMS132)   ---#
#---   Class 2: statistical inference              ---#
#-----------------------------------------------------#

#-- Example: Daily temperatures.

#-- data:
data=c(1.11, 2.87, 6.16, 2.48, 1.32, 4.26, 0.04, 2.00)

#-- some description od the data:
n=length(data)  # number of observations
summary(data)   # some summary of the data
hist(data, breaks=8, freq=FALSE, ylim=c(0,0.5), xlim=c(-1,10))

#-- a)

#-- b): unknown parameters!
pnorm(1.6, mean=mu, sd=sqrt(sigma2/n))

#-- c): variance=4
#   option 1: we can try to guess the distribution that generated our data:
plot(data, rep(0,n), xlim=c(-10,20), ylim=c(0,0.3),  pch=20, 
     main="data distribution", ylab="")
curve(dnorm(x, mean=0, sd=sqrt(4)),  col="red", lwd=3, add=TRUE)
curve(dnorm(x, mean=5, sd=sqrt(4)),  col="blue", lwd=3, add=TRUE)

#   option 2: the mean is interpreted as the averaged temperatures in several days.
mu=mean(data)
curve(dnorm(x, mean=mu, sd=sqrt(4)),  add=TRUE, col="black", lwd=3)


#-- d) distribution for mu
curve(dnorm(x, mean=0, sd=sqrt(100)), from=-20, to=20, ylim=c(0,0.1), lwd=3,
      xlab=expression(theta), ylab="", main="belief")


#-- e) we use Bayes theorem!
curve(dnorm(x, mean=0, sd=sqrt(100)), from=-20, to=20, ylim=c(0,0.2), lwd=3,
      xlab=expression(theta), ylab="", main="belief and updated belief") # our belief
updatedmean=( n*mean(data)/4 )/( n/4+1/100 )
updatedvariance=1/(n/4+1/100)
curve(dnorm(x, mean=updatedmean, sd=sqrt(updatedvariance)), from=-20, 
      to=20, ylim=c(0,0.1), lwd=3, add=TRUE, col="red") # updated belief

#-- f)
pnorm(3.21,mean=0, sd=sqrt(100))- pnorm(1.26,mean=0, sd=sqrt(100)) # our belief
pnorm(3.21,mean=updatedmean, sd=sqrt(updatedvariance))- pnorm(1.26,mean=updatedmean, sd=sqrt(updatedvariance)) # updated belief