Newer
Older
## 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
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)