Delete Conjugate_priors_Bayes.R

Conjugate_priors_Bayes.R

-## 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)
\ No newline at end of file