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

Update Conjugate_prior_and_RSTAN_MCMC.R

parent 1b86a87e
No related branches found
No related tags found
No related merge requests found
......@@ -36,78 +36,3 @@ conj <- bayesLMConjugate(y~x,data,n.samples,beta.prior.mean,beta.prior.precision
print(c(mean(conj$p.beta.tauSq.samples[,1]),sd((conj$p.beta.tauSq.samples[,1]))))
print(c(mean(conj$p.beta.tauSq.samples[,2]),sd((conj$p.beta.tauSq.samples[,2]))))
print(c(mean(sqrt(conj$p.beta.tauSq.samples[,3])),sd((sqrt(conj$p.beta.tauSq.samples[,3])))))
## For clarity, MCMC approach in R. Require RSTAN
## "0.7071068" corresponds to sqrt(0.5) (square root of the variance 0.5). Change it if prior
## precision matrix is changed. "b" in prior sigma is converted to rate rather than scale (1/scale)
## to match the conjugate approach package.
library(ggplot2)
library(rstan)
mymodel <- "
data {
int n;
vector[n] y;
vector[n] x;
}
parameters{
real alpha;
real beta;
real<lower=0> sigma;
}
model {
sigma ~ inv_gamma(0.5,1/0.5);
beta ~ normal(0,sigma/0.7071068);
alpha ~ normal(1,sigma/0.7071068);
y ~ normal(beta + x * alpha, sigma);
}
"
## Algorithm is NUTS as in PyMC3.
results <- stan(model_code = mymodel, data=list(n=14,x=x,y=y),iter=10000,chain=4,algorithm="NUTS")
beta_stan <- extract(results,pars="beta")
alpha_stan <- extract(results,pars="alpha")
sigma_stan <- extract(results,pars="sigma")
beta <- beta_stan$beta
alpha <- alpha_stan$alpha
sigma <- sigma_stan$sigma
## Check if same outputs than conjugate approach
print(c(mean(beta),sd(beta)))
print(c(mean(alpha),sd(alpha)))
print(c(mean(sigma),sd(sigma)))
## Below is the Bayesian inference to compute posterior S
## 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,20000)
quantile(pri,probs=c(0.05,0.95))
model_T = alpha*pri + beta + rnorm(n=20000,0,sd = sigma)
## PlioMIP temp
t = 0.8
stdt = 1
#
#t = -2.2
#stdt = 0.4
gaussobs <- rnorm(n=800000,t,stdt)
quantile(gaussobs,probs=c(0.05,0.95))
weight <- dnorm(model_T,t,stdt)
weight <- weight/sum(weight)
posterior <- sample(size = 1000000,x = pri,prob=weight,replace=TRUE)
quantile(posterior,probs=c(0.05,0.95))
median(posterior)
hist(posterior)
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