testStan.jl 2.26 KB
 Fredrik Bagge Carlson committed Oct 13, 2015 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 ``````######### 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 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) ``````