Commit 9c47a071 authored by Fredrik Bagge Carlson's avatar Fredrik Bagge Carlson
Browse files

initial work on incorporating Stan.jl

parent 469572fa
######### Stan program example ###########
# module Tmp
const σw0 = 1.0
const σw = 1
const σv = 1.0
const theta0 = [0.5, 25, 8]
s2piσv = log(sqrt(2*pi) * σv)
ProjDir = pwd()
function f_sample(x::Vector, t::Int64)
c = 8*cos(1.2*(t-1))
@inbounds for i = 1:length(x)
x[i] = 0.5*x[i] + 25*x[i]./(1+x[i]^2) + c + σw*randn()
end
return x
end
f_sample(x::Float64, t::Int64) = 0.5*x + 25*x/(1+x^2) + 8*cos(1.2*(t-1)) + σw*randn()
T = 100
M = 1
x = Array(Float64,T)
y = Array(Float64,T)
x0 = 0
x[1] = σw*randn()
y[1] = σv*randn()
for t = 1:T-1
x[t+1] = f_sample(x[t],t)
y[t+1] = 0.05x[t+1]^2 + σv*randn()
end # t
using Stan, Mamba
odemodel ="
data {
int<lower=1> T;
real y[T];
}
parameters {
real theta[3];
real sigma[2];
real x[T];
}
model {
sigma[1] ~ cauchy(1,1);
sigma[2] ~ cauchy(1,1);
theta[1] ~ cauchy(0.5,0.5);
theta[2] ~ cauchy(25,25);
theta[3] ~ cauchy(8,8);
x[1] ~ normal(0,1);
y[1] ~ normal(0.05*x[1]*x[1], sigma);
for (t in 1:(T-1)){
x[t+1] ~ normal(theta[1]*x[t] + theta[2]*x[t]/(1+x[t]*x[t]) + theta[3]*cos(1.2*(t-1)),sigma[1]);
y[t+1] ~ normal(0.05*x[t+1]*x[t+1], sigma[2]);
}
}
"
odedict =
Dict(
"T" => T,
"y" => y)
stanmodel = Stanmodel(name="ode", model=odemodel, nchains=4, update=10000)
@time sim1 = stan(stanmodel, [odedict], ProjDir, diagnostics=false, CmdStanDir=CMDSTAN_HOME)
## Subset Sampler Output to variables suitable for describe().
monitor = ["lp__", "accept_stat__"]
sim = sim1[1:size(sim1, 1), monitor, 1:size(sim1, 3)]
describe(sim)
println()
p = plot(sim, [:trace, :mean, :density, :autocor], legend=true)
draw(p, ncol=4, filename="$(stanmodel.name)-infoplot", fmt=:pdf)
## Subset Sampler Output to variables suitable for describe().
monitor = ["theta.1","theta.2","theta.3"]
sim = sim1[1:size(sim1, 1), monitor, 1:size(sim1, 3)]
describe(sim)
println()
p = plot(sim, [:trace, :mean, :density, :autocor], legend=true)
draw(p, ncol=4, filename="$(stanmodel.name)-thetaplot", fmt=:pdf)
## Subset Sampler Output to variables suitable for describe().
monitor = ["sigma.1", "sigma.2"]
sim = sim1[:, monitor, :]
describe(sim)
println()
p = plot(sim, [:trace, :mean, :density, :autocor], legend=true)
draw(p, nrow=4, ncol=4, filename="$(stanmodel.name)-sigmaplot", fmt=:pdf)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment