-
Notifications
You must be signed in to change notification settings - Fork 0
/
make_fig1.jl
executable file
·118 lines (105 loc) · 2.47 KB
/
make_fig1.jl
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
110
111
112
113
114
115
116
117
118
include("src/kalmanbucy.jl")
using .kalmanbucy
Ntrials = 100
T = 1000
dt = 0.001
par = [1., 2., 3.]
decpar = [10., sqrt(0.2), 3.]
rates = [0.03, 0.03, 0.0]
window = 20000 #window for moving averages
#----------------------------------------------------------------------
# Run simulation
#----------------------------------------------------------------------
t,X,dY,mu,P,af,bf,wf,mse = path(T,dt,par,decpar,rates)
#----------------------------------------------------------------------
# Produce figure
#----------------------------------------------------------------------
using Plots
fig1a = plot(;
xlabel = "time (in units of hidden state time-constant)",
grid = false,
ylimit = (-5, 5),
size = (1000, 300),
bottommargin = 1Plots.cm,
legend = nothing
)
plot!(
fig1a,
t[begin:100:end], X[begin:100:end];
color = :black, lw = 1,
label = "hidden state"
)
plot!(
fig1a,
t[begin:100:end], mu[begin:100:end];
color = :red, lw = 1,
label = "state estimate"
)
fig1b = plot(;
grid = false,
ylimit = (-5, 5),
size = (500, 300),
bottommargin = 1Plots.cm,
legend = nothing
)
plot!(
fig1b,
t[begin:100:10000], X[begin:100:10000];
color = :black, lw = 1,
label = "hidden state"
)
plot!(
fig1b,
t[begin:100:10000], mu[begin:100:10000];
color = :red, lw = 1,
label = "state estimate"
)
plot!(
fig1b,
t[begin:100:10000],
mu[begin:100:10000] - sqrt.(P[begin:100:10000]);
fill = mu[begin:100:10000] + sqrt.(P[begin:100:10000]),
color = :red, lw = 1, opacity = 0.3,
label = nothing
)
fig1c = plot(;
grid = false,
ylimit = (-5, 5),
size = (500, 300),
bottommargin = 1Plots.cm,
legend = nothing
)
plot!(
fig1c,
t[end-10000:100:end], X[end-10000:100:end];
color = :black, lw = 1,
label = "hidden state"
)
plot!(
fig1c,
t[end-10000:100:end], mu[end-10000:100:end];
color = :red, lw = 1,
label = "state estimate"
)
plot!(
fig1c,
t[end-10000:100:end],
mu[end-10000:100:end] - sqrt.(P[end-10000:100:end]);
fill = mu[end-10000:100:end] + sqrt.(P[end-10000:100:end]),
color = :red, lw = 1, opacity = 0.3,
label = nothing
)
fig1l = plot(; axis = false, grid = false, size = (200, 300))
plot!(fig1l, []; lw = 2, color = :black, label = " hidden state")
plot!(fig1l, []; lw = 2, color = :red, label = " state estimate")
fig1bottom = plot(
fig1b, fig1l, fig1c;
size = (1000, 200),
layout = grid(1, 3, widths = [0.45, 0.1, 0.45])
)
fig1 = plot(
fig1a, fig1bottom;
layout = (2, 1),
size = (1000, 500)
)
savefig(fig1, "fig/fig1.png")