#######################################################
###### KEY FUNCTIONS FROM THE "NCLBAYES" PACKAGE ######
#######################################################

##########################################################################
## If you're having problems installing the nclbayes package, fear not! ##
## Simply copy-and-paste the contents of this file into your RStudio    ##
## space and you can by-pass the package altogether!                    ##
##########################################################################

#####################################################################################
#### 1.  Highest density intervals for the Gamma, Beta and Inv-Chi distributions ####
#####################################################################################

### Firstly for a Gamma distribution:
hdiGamma = function(p = 0.95, a, b) 
{
    if (a <= 1) {
        lower = 0
        upper = qgamma(p, a, b)
    }
    else {
        zerofn = function(x) {
            (dgamma(qgamma(p + pgamma(x, a, b), a, b), a, b) - 
                dgamma(x, a, b))^2
        }
        maxl = qgamma(1 - p, a, b)
        res = optimize(zerofn, interval = c(0, maxl))
        lower = res$minimum
        upper = qgamma(p + pgamma(lower, a, b), a, b)
    }
    c(lower, upper)
}

### Now for a Beta distribution:
hdiBeta = function(p = 0.95, a, b) 
{
    if (a == b) {
        lower = qbeta((1 - p)/2, a, b)
        upper = qbeta((1 + p)/2, a, b)
    }
    else if ((a <= 1) & (b > 1)) {
        lower = 0
        upper = qbeta(p, a, b)
    }
    else if ((a > 1) & (b <= 1)) {
        lower = qbeta(1 - p, a, b)
        upper = 1
    }
    else {
        zerofn = function(x) {
            (dbeta(qbeta(p + pbeta(x, a, b), a, b), a, b) - dbeta(x, 
                a, b))^2
        }
        maxl = qbeta(1 - p, a, b)
        res = optimize(zerofn, interval = c(0, maxl))
        lower = res$minimum
        upper = qbeta(p + pbeta(lower, a, b), a, b)
    }
    c(lower, upper)
}

### Now for an Inv-chi distribution:
hdiInvchi = function(p = 0.95, a, b) 
{
    zerofn = function(x) {
        (dinvchi(qinvchi(p + pinvchi(x, a, b), a, b), a, b) - 
            dinvchi(x, a, b))^2
    }
    maxl = qinvchi(1 - p, a, b)
    res = optimize(zerofn, interval = c(0, maxl))
    lower = res$minimum
    upper = qinvchi(p + pinvchi(lower, a, b), a, b)
    c(lower, upper)
}

#########################################################
#### 2.  Functions for the Normal-Gamma distribution ####
#########################################################

### Firstly, a function to compute the joint density of a N-Ga:
dnormgamma = function(mu, tau, b, c, g, h) 
{
    dnorm(mu, b, 1/sqrt(c * tau)) * dgamma(tau, g, h)
}

### Now, the function to produce contour plots for the N-Ga:
 NGacontour = function(mu, tau, b, c, g, h, p = NULL, ...) 
{
    pdf = function(mu, tau) {
        dnormgamma(mu, tau, b, c, g, h)
    }
    z = outer(mu, tau, pdf)
    if (length(p) == 0) {
        contour(mu, tau, z, xlab = expression(mu), ylab = expression(tau), 
            ...)
    }
    else {
        NGaclevels = function(p, b, c, g, h) {
            n = 1e+05
            tau = rgamma(n, g, h)
            mu = rnorm(n, b, 1/sqrt(c * tau))
            y = dnormgamma(mu, tau, b, c, g, h)
            y = sort(y)
            signif(y[floor(n * (1 - p))], 2)
        }
        pdflevels = NGaclevels(p, b, c, g, h)
        contour(mu, tau, z, levels = pdflevels, xlab = expression(mu), 
            ylab = expression(tau), ...)
    }
}


######################################################################################################
#### 3.  Functions to compute densities, probabilities and quantiles for some "new" distributions ####
######################################################################################################

### The generalised t-distribution: density, probabilities and quantiles:
dgt = function(x, a, b, c) 
{
    dt((x - b)/sqrt(c), a)/sqrt(c)
}

pgt = function (x, a, b, c) 
{
    pt((x - b)/sqrt(c), a)
}

qgt = function(p, a, b, c) 
{
    b + sqrt(c) * qt(p, a)
}

### The inverse-Chi distribution: density, probabilities and quantiles:
dinvchi = function(x, a, b) 
{
    2 * dgamma(1/x^2, a, b)/x^3
}

pinvchi = function(x, a, b) 
{
    1 - pgamma(1/x^2, a, b)
}

qinvchi = function(p, a, b) 
{
    1/sqrt(qgamma(1 - p, a, b))
}


###################################
#### 4. MCMC functions ############
###################################

##Lee's code to perform a Gibbs Sampler for Example 4.2 i the notes:

gibbs2 = function(N, data.values=c(100,15,4.5), prior.params=c(10,1/100, 3,12)){
  n = data.values[1]
  xbar = data.values[2]
  s = data.values[3]
  b = prior.params[1]
  c = prior.params[2]
  g = prior.params[3]
  h = prior.params[4]
  mu = vector("numeric", N)
  tau = vector("numeric", N)
  mu[1] = b
  tau[1] = g/h
  for(j in 2:N){
    print(j)
    mu[j] = rnorm(1, (b*c+n*tau[j-1]*xbar)/(c+n*tau[j-1]), sqrt(1/(c+n*tau[j-1])))
    tau[j] = rgamma(1, g+(n/2), h+(n/2)*(s^2+((xbar-mu[j])^2)))}
  results = matrix(ncol=2, nrow=N)
  results[,1] = mu
  results[,2] = tau
  return(results)}

##The equivalent code from the "nclbayes" package:

gibbsNormal = function(N, initial, priorparam, n, xbar, s) 
{
    b = priorparam[1]
    c = priorparam[2]
    g = priorparam[3]
    h = priorparam[4]
    output = matrix(ncol = 2, nrow = N)
    colnames(output) = c("mu", "tau")
    mu = initial[1]
    tau = initial[2]
    output[1, ] = c(mu, tau)
    for (i in 2:N) {
        muprec = c + n * tau
        mumean = (b * c + n * tau * xbar)/muprec
        mu = rnorm(1, mumean, 1/sqrt(muprec))
        tauh = h + 0.5 * n * (s^2 + (xbar - mu)^2)
        tau = rgamma(1, g + n/2, tauh)
        output[i, ] = c(mu, tau)
    }
    output
}

#Metropolis-within-Gibbs sampling:
mwgGamma =function(N, initial, innov, priorparam, n, xbar, xgbar, show = TRUE) 
{
    output = matrix(ncol = 2, nrow = N)
    colnames(output) = c("alpha", "lambda")
    a = priorparam[1]
    b = priorparam[2]
    c = priorparam[3]
    d = priorparam[4]
    alpha = initial
    naccept = 0
    for (i in 1:N) {
        lambda = rgamma(1, c + n * alpha, d + n * xbar)
        proposal = rnorm(1, alpha, innov)
        if (proposal > 0) {
            logaprob = (a - 1) * log(proposal/alpha) + n * log(gamma(alpha)/gamma(proposal)) + 
                (b + n * log(xgbar) + n * log(lambda)) * (proposal - 
                  alpha)
            u = runif(1)
            if (log(u) < logaprob) {
                alpha = proposal
                naccept = naccept + 1
            }
        }
        output[i, ] = c(alpha, lambda)
    }
    if (show == TRUE) 
        message(paste("Acceptance rate = ", naccept/N))
    output
}

#Processing MCMC output:
mcmcProcess = function(input, burnin = 1000, thin = 1) 
{
    output = input[seq(burnin + 1, nrow(input), by = thin), ]
    if (ncol(input) == 1) {
        as.matrix(output)
    }
    else {
        output
    }
}

mcmcAnalysis = function(mat, rows = 4, lag.max = 100, bins = 30, show = TRUE, 
    plot = TRUE) 
{
    d = dim(mat)
    p = d[2]
    summ = summary(mat)
    if (show == TRUE) {
        message(paste("N =", d[1], "iterations"))
        print(summ)
        message("Standard deviations:")
        print(apply(mat, 2, sd))
    }
    if (plot == TRUE) {
        names = colnames(mat)
        op = par(mfrow = c(rows, 3))
        for (i in 1:p) {
            plot(ts(mat[, i]), main = names[i], ylab = "Value", 
                xlab = "Iteration")
            acf(mat[, i], lag.max = lag.max, main = names[i], 
                ci = 0)
            hist(mat[, i], bins, main = names[i], xlab = "Value", 
                freq = FALSE)
        }
        par(op)
    }
    invisible(summ)
}
mcmcCi = function(input, level = 0.95) 
{
    N = nrow(input)
    p = ncol(input)
    output = matrix(ncol = 2, nrow = p)
    rownames(output) = colnames(input)
    for (i in 1:p) {
        x = sort(input[, i])
        output[i, 1] = x[floor(N * (1 - level)/2)]
        output[i, 2] = x[floor(N * (1 + level)/2)]
    }
    output
}