testStan.jl 2.26 KB
Newer Older
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<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)