# R-code for obtaining ML estimate of sigma by using quadratic regression
# Output from ESM S1-1 is used.

# Create the data frame of histogram
Strans <- out$sims.list$Strans
freq <- hist(Strans, breaks=70, plot=FALSE)
histlength <- length(freq$breaks)
x <- (freq$breaks[1:(histlength-1)] + freq$breaks[2:histlength])/2
histStrans <- data.frame(x, y=freq$counts)

# Extract the central part of histogram
histStrans <- subset(histStrans, y > mean(y))
histStrans$logy <- log(histStrans$y)

# Quadratic regression
histStrans.glm <- glm(logy ~ x + I(x^2), data=histStrans)

# Estimation of mode
estimate <- coef(histStrans.glm)
a <- estimate["(Intercept)"]
b <- estimate["x"]
c <- estimate["I(x^2)"]
StransMode <- -b/(2*c)

# back-transformaion of sigma
lambdaS <- -0.5
sigma <- (lambdaS*StransMode+1)^(1/lambdaS)
cat("ML estimate by quadratic regression =", sigma, "\n")

# Plot the log-likelihood and the quadratic curve
plot(logy~x, data=histStrans, xlab="Transformed sigma", ylab="Logarithmic likelihood", cex.lab=1.2, cex.axis=1)
quad <- function(x) a + b*x + c*x^2
plot(quad, -6.5, -2.5, add=TRUE)

# Emprical Jeffreys prior estimate (EJP)
abline(v=-4.61100, lty=1)
text(-4.48, 7.2, "EJP")

# dclone estimate (DC)
dclone.sigma <- (0.08122681^lambdaS-1)/lambdaS
abline(v=dclone.sigma, lty=2)
text(-5.12, 7.2, "DC")