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

Add new file

parent 8ba22aa8
No related branches found
No related tags found
No related merge requests found
## 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.
## Priors on the parameters can be freely changed just as when using PyMC3. Please refer to
## RSTAN for a list of possible prior distributions.
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)
\ 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