Skip to content
Snippets Groups Projects
Commit 89d2ef8d authored by Martin Renoult's avatar Martin Renoult
Browse files

Add new file

parent 70df8d03
No related branches found
No related tags found
No related merge requests found
## 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 <- 0.5
prior.rate <- 0.5
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])
sigma <- as.vector(bmod$p.beta.tauSq.samples[,3])
hist(beta)
hist(alpha)
hist(sigma)
write.csv(beta, file="intercept_conj.csv")
write.csv(alpha, file="slope_conj.csv")
write.csv(sigma,file="sigma_conj.csv")
## 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)
\ No newline at end of file
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment