Skip to content

Commit 8cf5126

Browse files
compopscompops
committed
refactorizated R code and re-run it for new paper draft.
Former-commit-id: acf4712
1 parent e8890bc commit 8cf5126

24 files changed

Lines changed: 1221 additions & 626 deletions

‎.vscode/.ropeproject/config.py‎

Lines changed: 100 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,100 @@
1+
# The default ``config.py``
2+
# flake8: noqa
3+
4+
5+
def set_prefs(prefs):
6+
"""This function is called before opening the project"""
7+
8+
# Specify which files and folders to ignore in the project.
9+
# Changes to ignored resources are not added to the history and
10+
# VCSs. Also they are not returned in `Project.get_files()`.
11+
# Note that ``?`` and ``*`` match all characters but slashes.
12+
# '*.pyc': matches 'test.pyc' and 'pkg/test.pyc'
13+
# 'mod*.pyc': matches 'test/mod1.pyc' but not 'mod/1.pyc'
14+
# '.svn': matches 'pkg/.svn' and all of its children
15+
# 'build/*.o': matches 'build/lib.o' but not 'build/sub/lib.o'
16+
# 'build//*.o': matches 'build/lib.o' and 'build/sub/lib.o'
17+
prefs['ignored_resources'] = ['*.pyc', '*~', '.ropeproject',
18+
'.hg', '.svn', '_svn', '.git', '.tox']
19+
20+
# Specifies which files should be considered python files. It is
21+
# useful when you have scripts inside your project. Only files
22+
# ending with ``.py`` are considered to be python files by
23+
# default.
24+
#prefs['python_files'] = ['*.py']
25+
26+
# Custom source folders: By default rope searches the project
27+
# for finding source folders (folders that should be searched
28+
# for finding modules). You can add paths to that list. Note
29+
# that rope guesses project source folders correctly most of the
30+
# time; use this if you have any problems.
31+
# The folders should be relative to project root and use '/' for
32+
# separating folders regardless of the platform rope is running on.
33+
# 'src/my_source_folder' for instance.
34+
#prefs.add('source_folders', 'src')
35+
36+
# You can extend python path for looking up modules
37+
#prefs.add('python_path', '~/python/')
38+
39+
# Should rope save object information or not.
40+
prefs['save_objectdb'] = True
41+
prefs['compress_objectdb'] = False
42+
43+
# If `True`, rope analyzes each module when it is being saved.
44+
prefs['automatic_soa'] = True
45+
# The depth of calls to follow in static object analysis
46+
prefs['soa_followed_calls'] = 0
47+
48+
# If `False` when running modules or unit tests "dynamic object
49+
# analysis" is turned off. This makes them much faster.
50+
prefs['perform_doa'] = True
51+
52+
# Rope can check the validity of its object DB when running.
53+
prefs['validate_objectdb'] = True
54+
55+
# How many undos to hold?
56+
prefs['max_history_items'] = 32
57+
58+
# Shows whether to save history across sessions.
59+
prefs['save_history'] = True
60+
prefs['compress_history'] = False
61+
62+
# Set the number spaces used for indenting. According to
63+
# :PEP:`8`, it is best to use 4 spaces. Since most of rope's
64+
# unit-tests use 4 spaces it is more reliable, too.
65+
prefs['indent_size'] = 4
66+
67+
# Builtin and c-extension modules that are allowed to be imported
68+
# and inspected by rope.
69+
prefs['extension_modules'] = []
70+
71+
# Add all standard c-extensions to extension_modules list.
72+
prefs['import_dynload_stdmods'] = True
73+
74+
# If `True` modules with syntax errors are considered to be empty.
75+
# The default value is `False`; When `False` syntax errors raise
76+
# `rope.base.exceptions.ModuleSyntaxError` exception.
77+
prefs['ignore_syntax_errors'] = False
78+
79+
# If `True`, rope ignores unresolvable imports. Otherwise, they
80+
# appear in the importing namespace.
81+
prefs['ignore_bad_imports'] = False
82+
83+
# If `True`, rope will insert new module imports as
84+
# `from <package> import <module>` by default.
85+
prefs['prefer_module_from_imports'] = False
86+
87+
# If `True`, rope will transform a comma list of imports into
88+
# multiple separate import statements when organizing
89+
# imports.
90+
prefs['split_imports'] = False
91+
92+
# If `True`, rope will sort imports alphabetically by module name
93+
# instead of alphabetically by import statement, with from imports
94+
# after normal imports.
95+
prefs['sort_imports_alphabetically'] = False
96+
97+
98+
def project_opened(project):
99+
"""This function is called after opening the project"""
100+
# Do whatever you like here!

‎.vscode/.ropeproject/objectdb‎

6 Bytes
Binary file not shown.

‎r/example1-lgss-plotData.R‎

Lines changed: 143 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,143 @@
1+
##############################################################################
2+
#
3+
# Example of fully-adapted particle filtering
4+
# in a linear Gaussian state space model
5+
#
6+
#
7+
# Copyright (C) 2017 Johan Dahlin < liu (at) johandahlin.se >
8+
#
9+
# This program is free software; you can redistribute it and/or modify
10+
# it under the terms of the GNU General Public License as published by
11+
# the Free Software Foundation; either version 2 of the License, or
12+
# (at your option) any later version.
13+
#
14+
# This program is distributed in the hope that it will be useful,
15+
# but WITHOUT ANY WARRANTY; without even the implied warranty of
16+
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17+
# GNU General Public License for more details.
18+
#
19+
# You should have received a copy of the GNU General Public License along
20+
# with this program; if not, write to the Free Software Foundation, Inc.,
21+
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
22+
#
23+
##############################################################################
24+
25+
# Import helper
26+
source("stateEstimationHelper.R")
27+
28+
# Set the random seed to replicate results in tutorial
29+
set.seed(10)
30+
31+
# Save plot to file
32+
savePlotToFile <- TRUE
33+
34+
##############################################################################
35+
# Define the model
36+
##############################################################################
37+
38+
# Here, we use the following model
39+
#
40+
# x[tt+1] = phi * x[tt] + sigmav * v[tt]
41+
# y[tt] = x[tt] + sigmae * e[tt]
42+
#
43+
# where v[tt] ~ N(0,1) and e[tt] ~ N(0,1)
44+
45+
# Set the parameters of the model
46+
phi <- 0.75
47+
sigmav <- 1.00
48+
sigmae <- 0.10
49+
50+
# Set the number of time steps to simulate
51+
T <- 250
52+
53+
# Set the initial state
54+
initialState <- 0
55+
56+
##############################################################################
57+
# Generate data
58+
##############################################################################
59+
60+
data <- generateData(c(phi, sigmav, sigmae), T, initialState)
61+
x <- data$x
62+
y <- data$y
63+
64+
##############################################################################
65+
# Plotting
66+
##############################################################################
67+
68+
# Export plot to file
69+
if (savePlotToFile) {
70+
cairo_pdf("figures/lgss-data.pdf",
71+
height = 3,
72+
width = 8)
73+
}
74+
75+
grid = seq(0, T)
76+
77+
# Plot the latent state and observations
78+
layout(matrix(1:3, 1, 3, byrow = TRUE))
79+
par(mar = c(4, 5, 0, 0))
80+
81+
plot(
82+
grid,
83+
x,
84+
col = "#D95F02",
85+
lwd = 1,
86+
type = "l",
87+
xlab = "time",
88+
ylab = expression("latent state " * x[t]),
89+
bty = "n",
90+
ylim = c(-4, 6)
91+
)
92+
polygon(c(grid, rev(grid)),
93+
c(x, rep(-4, T + 1)),
94+
border = NA,
95+
col = rgb(t(col2rgb("#D95F02")) / 256, alpha = 0.25))
96+
97+
plot(
98+
grid[-1],
99+
y[-1],
100+
col = "#1B9E77",
101+
lwd = 1,
102+
type = "l",
103+
xlab = "time",
104+
ylab = expression("observation " * y[t]),
105+
bty = "n",
106+
ylim = c(-4, 6)
107+
)
108+
polygon(c(grid[-1], rev(grid[-1])),
109+
c(y[-1], rep(-4, T)),
110+
border = NA,
111+
col = rgb(t(col2rgb("#1B9E77")) / 256, alpha = 0.25))
112+
113+
foo = acf(y[-1], plot = F, lag.max = 25)
114+
115+
plot(
116+
foo$lag,
117+
foo$acf,
118+
col = "#66A61E",
119+
lwd = 1.5,
120+
type = "l",
121+
xlab = "time",
122+
ylab = expression("ACF of " * y[t]),
123+
bty = "n",
124+
ylim = c(-0.2, 1),
125+
xlim = c(0, 25)
126+
)
127+
polygon(
128+
c(foo$lag, rev(foo$lag)),
129+
c(foo$acf, rep(0.0, length(foo$lag))),
130+
border = NA,
131+
col = rgb(t(col2rgb("#66A61E")) / 256, alpha = 0.25)
132+
)
133+
abline(h = -1.96 / sqrt(T), lty = "dotted")
134+
abline(h = 1.96 / sqrt(T), lty = "dotted")
135+
136+
137+
if (savePlotToFile) {
138+
dev.off()
139+
}
140+
141+
###################################################################################
142+
# End of file
143+
###################################################################################

‎r/example1-lgss.R‎

Lines changed: 69 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -28,6 +28,12 @@ source("stateEstimationHelper.R")
2828
# Set the random seed to replicate results in tutorial
2929
set.seed(10)
3030

31+
# Should the results be loaded from file (to quickly generate plots)
32+
loadSavedWorkspace <- FALSE
33+
34+
# Save plot to file
35+
savePlotToFile <- TRUE
36+
3137

3238
##############################################################################
3339
# Define the model
@@ -51,6 +57,7 @@ T <- 250
5157
# Set the initial state
5258
initialState <- 0
5359

60+
5461
##############################################################################
5562
# Generate data
5663
##############################################################################
@@ -59,6 +66,13 @@ data <- generateData(c(phi, sigmav, sigmae), T, initialState)
5966
x <- data$x
6067
y <- data$y
6168

69+
# Export plot to file
70+
if (savePlotToFile) {
71+
cairo_pdf("figures/example1-lgss.pdf",
72+
height = 10,
73+
width = 8)
74+
}
75+
6276
# Plot the latent state and observations
6377
layout(matrix(c(1, 1, 2, 2, 3, 4), 3, 2, byrow = TRUE))
6478
par (mar = c(4, 5, 0, 0))
@@ -75,16 +89,28 @@ plot(
7589
ylim = c(-6, 6),
7690
bty = "n"
7791
)
92+
polygon(c(grid, rev(grid)),
93+
c(y, rep(-6, T + 1)),
94+
border = NA,
95+
col = rgb(t(col2rgb("#1B9E77")) / 256, alpha = 0.25))
96+
7897

7998
###################################################################################
8099
# State estimation using the particle filter and Kalman filter
81100
###################################################################################
82101

83-
# Using noParticles = 20 particles and plot the estimate of the latent state
84-
noParticles <- 20
85-
outputPF <- particleFilter(y, c(phi, sigmav, sigmae), noParticles, initialState)
86-
outputKF <- kalmanFilter(y, c(phi, sigmav, sigmae), initialState, 0.01)
87-
difference <- outputPF$xHatFiltered - outputKF$xHatFiltered[-(T + 1)]
102+
if (loadSavedWorkspace) {
103+
load("savedWorkspaces/example1-lgss.RData")
104+
} else {
105+
# Using noParticles = 20 particles and plot the estimate of the latent state
106+
noParticles <- 20
107+
outputPF <-
108+
particleFilter(y, c(phi, sigmav, sigmae), noParticles, initialState)
109+
outputKF <-
110+
kalmanFilter(y, c(phi, sigmav, sigmae), initialState, 0.01)
111+
difference <-
112+
outputPF$xHatFiltered - outputKF$xHatFiltered[-(T + 1)]
113+
}
88114

89115
grid <- seq(0, T - 1)
90116
plot(
@@ -97,18 +123,25 @@ plot(
97123
ylim = c(-0.1, 0.1),
98124
bty = "n"
99125
)
126+
polygon(
127+
c(grid, rev(grid)),
128+
c(difference, rep(-0.1, T)),
129+
border = NA,
130+
col = rgb(t(col2rgb("#7570B3")) / 256, alpha = 0.25)
131+
)
100132

101133
# Compute bias and MSE
102134
logBiasMSE <- matrix(0, nrow = 7, ncol = 2)
103135
gridN <- c(10, 20, 50, 100, 200, 500, 1000)
104136

105137
for (ii in 1:length(gridN)) {
106-
pfEstimate <- particleFilter(y, c(phi, sigmav, sigmae), gridN[ii], initialState)
138+
pfEstimate <-
139+
particleFilter(y, c(phi, sigmav, sigmae), gridN[ii], initialState)
107140
pfEstimate <- pfEstimate$xHatFiltered
108141
kfEstimate <- outputKF$xHatFiltered[-(T + 1)]
109142

110-
logBiasMSE[ii, 1] <-log(mean(abs(pfEstimate - kfEstimate)))
111-
logBiasMSE[ii, 2] <-log(mean((pfEstimate - kfEstimate)^2))
143+
logBiasMSE[ii, 1] <- log(mean(abs(pfEstimate - kfEstimate)))
144+
logBiasMSE[ii, 2] <- log(mean((pfEstimate - kfEstimate) ^ 2))
112145
}
113146

114147
# Plot the bias and MSE for comparison
@@ -119,9 +152,15 @@ plot(
119152
type = "l",
120153
xlab = "no. particles (N)",
121154
ylab = "log-bias",
122-
ylim = c(-7, -3),
155+
ylim = c(-7,-3),
123156
bty = "n"
124157
)
158+
polygon(
159+
c(gridN, rev(gridN)),
160+
c(logBiasMSE[, 1], rep(-7, length(gridN))),
161+
border = NA,
162+
col = rgb(t(col2rgb("#E7298A")) / 256, alpha = 0.25)
163+
)
125164
points(gridN, logBiasMSE[, 1], col = "#E7298A", pch = 19)
126165

127166
plot(
@@ -132,13 +171,32 @@ plot(
132171
type = "l",
133172
xlab = "no. particles (N)",
134173
ylab = "log-MSE",
135-
ylim = c(-12, -6),
174+
ylim = c(-12,-6),
136175
bty = "n"
137176
)
177+
polygon(
178+
c(gridN, rev(gridN)),
179+
c(logBiasMSE[, 2], rep(-12, length(gridN))),
180+
border = NA,
181+
col = rgb(t(col2rgb("#66A61E")) / 256, alpha = 0.25)
182+
)
138183
points(gridN, logBiasMSE[, 2], col = "#66A61E", pch = 19)
139184

185+
# Close the plotting device
186+
if (savePlotToFile) {
187+
dev.off()
188+
}
189+
190+
191+
##############################################################################
192+
# Compute and save the results
193+
##############################################################################
194+
140195
# Save the workspace to file
141-
save("example1-lgss.RData")
196+
if (!loadSavedWorkspace) {
197+
save.image("savedWorkspaces/example1-lgss.RData")
198+
}
199+
142200

143201
###################################################################################
144202
# End of file

0 commit comments

Comments
 (0)