source("hdr.txt")
priormean <- 300
priorsd <- 160
nu0 <- 2*priormean^2/priorsd^2+4
nu0 <- round(nu0)
S0 <- priormean*(nu0-2)
theta0 <- 110
n0 <- 15
x=c(141,102,73,171,137,91,81,157,146,69,121,134)
n <- length(x)
xbar <- mean(x)
S <- sum((x-xbar)^2)
s <- sd(x)
nu1 <- nu0+n
n1 <- n0+n
theta1 <- (n0*theta0+n*xbar)/n1
S1 <- S0 + S + (n0^(-1)+n^(-1))^(-1)*(theta0-xbar)^2
s <- sqrt(S1/nu1)
tvalues <- ht(0.95,nu1)
cat("Interval for mu                                   ",
    (s/sqrt(n1))*tvalues+theta1,"\n")
cat("Posterior mean of phi                             ",
    S1/(nu1-2),"\n")
cat("Values of chi^2 corresponding to HDR for log chi^2",
    hchisq(0.95,nu1),"\n")
cat("Interval for phi                                  ",
    sort(S1/hchisq(0.95,19)),"\n")