## Library: In particular, spBayes is used for the conjugate problem

library(coda)
library(magic)
library(abind)
library(Formula)
library(Matrix)
library(spBayes)
library(extraDistr)

## PlioMIP1 + 2 ensemble
x <- c(3.2, 3.4, 4.05, 2.8, 4.1, 3.2, 3.3, 3.1,
       3.37, 2.6,4.50, 2.29,5.3,4.3)

y <- c(1.0256042, 1.3337059, 1.9900227, 1.1576538, 2.1806774, 1.151104, 1.9331722, 1.445343,
       2.14, 0.9211941,2.1174774, 1.3736095,3.4950447,2.94250999999997)

model_data <- data.frame(y,x)

## p defines number of parameters (2: Intercept and Slope)
n <- nrow(model_data)
p <- 2
n.samples <- 10000

beta.prior.mean <- c(0,1)

beta.prior.precision <- 1*diag(p)

## Emulate the original Cauchy prior with some specific parameters of InvGamma
prior.shape <- 2
prior.rate <- 2

bmod <- bayesLMConjugate(y ~ x, data=model_data, n.samples,beta.prior.mean, beta.prior.precision, prior.shape, prior.rate)
beta <- as.vector(bmod$p.beta.tauSq.samples[,1])
alpha <- as.vector(bmod$p.beta.tauSq.samples[,2])
tauSq <- as.vector(bmod$p.beta.tauSq.samples[,3])
sigma <- sqrt(tauSq)

hist(beta)
hist(alpha)
hist(sigma)

## Truncated Cauchy prior
trunc <- function(loc, scale, a, b, sz) {
ua = atan((a - loc)/scale)/pi + 0.5
ub = atan((b - loc)/scale)/pi + 0.5
U <- runif(n = sz,ua,ub)
rvs <<- loc + scale * tan(pi*(U-0.5))}
pri <- trunc(2.5,3,1/Inf,Inf,100000)

quantile(pri,probs=c(0.05,0.95))

model_T = alpha*pri + beta + rnorm(n=100000,0,sd = sigma)

## PlioMIP temp
t = 0.8
stdt = 1

gaussobs <- rnorm(n=8000000,t,stdt)
quantile(gaussobs,probs=c(0.05,0.95))

weight <- dnorm(model_T,t,stdt)
weight <- weight/sum(weight)

posterior <- sample(size = 10000,x = pri,prob=weight)
quantile(posterior,probs=c(0.05,0.95))
median(posterior)
hist(posterior)