# Right/Wrong assumptions and Bayesian intervals.

# a)
x <-c(1.10, 0.87, 0.84, 2.62, 17.37, 7.38, 3.24, 5.74, 0.88, 8.34, 4.57,
      7.43, 26.54, 6.33, 6.21, 11.26, 3.93, 2.02, 3.53, 14.19, 3.85,
      1.76, 3.40, 0.64, 0.36)

n <- length(x)


hist(x)


# b) normal assumption:
gamma=0.9
aintN <- mean(x) - qt((1+gamma)/2, n-1)*sd(x)/sqrt(n)
bintN <- mean(x) + qt((1+gamma)/2, n-1)*sd(x)/sqrt(n)
c(aintN, bintN)

bintN-aintN

# lower bounded interval:
mean(x) - qt(gamma, n-1)*sd(x)/sqrt(n)

#  c) exponential assumption:
qchisq(0.05, 2*n)
qchisq(0.95, 2*n)
curve(dchisq(x, 2*n), 0, 100)

aintE <- 2*sum(x)/qchisq(0.95, 2*n)
bintE <- 2*sum(x)/qchisq(0.05, 2*n)
c(aintE, bintE)

bintE-aintE


# Bayesian intervals:
# d)
alpha= 1
beta=0.25

curve(dchisq(x, 2*alpha), 0, 100)
qchisq(0.05, 2*alpha)
qchisq(0.95, 2*alpha)

2*beta/qchisq(0.95, 2*alpha)
2*beta/qchisq(0.05, 2*alpha)

curve(dgamma(x, 1, 0.1),0, 100)
curve(dgamma(x, 1, 0.1),0, 100, ylim=c(0,0.2))

# e)
curve(dchisq(x, 2*(alpha+n)), 0, 100)
qchisq(0.05, 2*(alpha+n))
qchisq(0.95, 2*(alpha+n))

2*(beta+sum(x))/qchisq(0.95, 2*(alpha+n))
2*(beta+sum(x))/qchisq(0.05, 2*(alpha+n))


# g)
alpha=1; beta=0.25
alpha/beta; alpha/beta^2
curve(dgamma(x, alpha, beta),0, 100, lwd=2)
# Prior region
c(2*beta/qchisq(0.95, 2*alpha), 2*beta/qchisq(0.05, 2*alpha))
# posterior credible region:
c(2*(beta+sum(x))/qchisq(0.95, 2*(alpha+n)), 2*(beta+sum(x))/qchisq(0.05, 2*(alpha+n)))

alpha=1; beta=1
curve(dgamma(x, alpha, beta),0, 100, add=TRUE, lwd=3, col='red')
# Prior region
c(2*beta/qchisq(0.95, 2*alpha), 2*beta/qchisq(0.05, 2*alpha))
# posterior credible region:
c(2*(beta+sum(x))/qchisq(0.95, 2*(alpha+n)), 2*(beta+sum(x))/qchisq(0.05, 2*(alpha+n)))

alpha=1; beta=0.01
curve(dgamma(x, alpha, beta),0, 100, add=TRUE, lwd=3, col='blue')
# Prior region
c(2*beta/qchisq(0.95, 2*alpha), 2*beta/qchisq(0.05, 2*alpha))
# posterior credible region:
c(2*(beta+sum(x))/qchisq(0.95, 2*(alpha+n)), 2*(beta+sum(x))/qchisq(0.05, 2*(alpha+n)))


# plot:
plot(0,0, xlim=c(0,20), ylim=c(0,4), type="n")
alpha=1; beta=0.25; alpha/beta
segments(2*beta/qchisq(0.95, 2*alpha), 1, 2*beta/qchisq(0.05, 2*alpha),1, lwd=3)
segments(2*(beta+sum(x))/qchisq(0.95, 2*(alpha+n)), 1.2, 2*(beta+sum(x))/qchisq(0.05, 2*(alpha+n)),1.2, lwd=3, col="red")

alpha=1; beta=1; alpha/beta # the mean is close to very small values
segments(2*beta/qchisq(0.95, 2*alpha), 2, 2*beta/qchisq(0.05, 2*alpha),2, lwd=3)
segments(2*(beta+sum(x))/qchisq(0.95, 2*(alpha+n)), 2.2, 2*(beta+sum(x))/qchisq(0.05, 2*(alpha+n)),2.2, lwd=3, col="red")

alpha=1; beta=0.01; alpha/beta
segments(2*beta/qchisq(0.95, 2*alpha), 3, 2*beta/qchisq(0.05, 2*alpha),3, lwd=3)
segments(2*(beta+sum(x))/qchisq(0.95, 2*(alpha+n)), 3.2, 2*(beta+sum(x))/qchisq(0.05, 2*(alpha+n)),3.2, lwd=3, col="red")

segments(2*sum(x)/qchisq(0.95, 2*n), 0.5, 2*sum(x)/qchisq(0.05, 2*n),0.5, lwd=3, col="green")
abline(v=c(2*sum(x)/qchisq(0.95, 2*n), 2*sum(x)/qchisq(0.05, 2*n)), lty=3)