[PYTHON] Try a state-space model (Jupyter Notebook + IR kernel)
Outline of state space model
The book is Iwanami Data Science Vol.1. ([Special feature] Bayesian inference and MCMC free software)
The conceptual diagram of the state space model is quoted from here.
A state-space model has a time-series ** "state" **, and based on this, ** "observed values" ** containing certain variations are generated. That is. Iwanami Data Science Vol.1 is a special feature of "Bayesian Inference and MCMC Free Software", but Chapter 3 is an article about state-space models (written by Mr. Ito), and other chapters also explain spatial models. come. Chapter 3 explained how to solve this with R and JAGS (Just Another Gibbs Sampler).
When I tried this, I initially tried it with Python + PyMC3, but it didn't work.
import pymc3
(Omitted)
model1 = pm.Model()
with model1: # model definition
alpha0 = pm.Uniform('alpha0', lower=0. , upper=10.)
sigma_y = pm.Uniform('sigma_y', lower=0., upper=100.)
sigma_alp = pm.Uniform('sigma_alp', lower=0., upper=100.)
for i in range(1,N): #The following two lines are ERROR.I wanted to define the state quantity α...
alpha[i] = pm.Normal(1, mu=alpha[i-1], sd=sigma[0])
y_obs = pm.Normal('y_obs', mu=alpha, sd=sigma[0], observed=ydat)
(Omitted)
I tried such a code, but in the for loop, I get an error when defining the state quantity model $ \ alpha_t = \ alpha_ {t-1} + \ eta $. If you think about it, ** pymc3.Normal () ** (** pm.Normal () ** in the list) is a class constructor for generating a model of normal distribution, so in its argument It is considered difficult to include the class (alpha [i-1]) to be generated. (I searched the Internet variously, but couldn't find a solution.)
I considered using ** PyStan ** and Kalman filter implementation ** PyKalman ** as other libraries that can handle such state-space models (or dynamic linear models), but this time I followed the book. I decided to use "R". The program used is as follows.
- R ver.3.2.2 --Jupyter Notebook (notebook is ver.4.0.6) + IRkernel ... This time, I tried IRkernel for the first time. --JAGS ver.4.0.0 + {rjags} package -{dlm} package
About Jupyter Notebook + IRkernel environment
Jupyter Notebook is an updated version of the IPython Notebook you know. By switching the language processing system and kernel, it is possible to handle other than Python. (It was named after Jupyter = Julia + Python + R.) This time, I installed and used IRkernel for R language.
The installation procedure is as follows.
-Install R language (if necessary). --Install IRkernel (https://github.com/IRkernel/IRkernel) under R environment.
install.packages(c('rzmq','repr','IRkernel','IRdisplay'),
repos = c('http://irkernel.github.io/', getOption('repos')))
IRkernel::installspec()
--Start with "> Jupyte Notebook" on the console. (At first, an error occurred due to lack of some R dependent packages, but it can be started by installing the R package as requested by the error message.)
** Fig. When Jupyter Notebook is started **
As shown above, you can select "R" from the new note creation menu.
Implementation of state space model by JAGS
Follow the "Modeling of Time Series / Spatial Data" chapter of the book "Iwanami Data Science Vol.1" to create the following notes.
As with previous IPython Notebooks, you can mix Markdown text with executable code. Also, from this update of Notebook, I'm happy to be able to render ** MathJax ** formulas offline.
One tip for plotting figures with R. (Notebook In [1]: Command described.)
library(repr)
# Change plot size to 4 x 3
options(repr.plot.width=4, repr.plot.height=3)
In the latter half of this note, we plot using {ggplot2}, but if nothing is specified, the figure will be generated to fit the full width of the Web Browser. The figure is quite large (although most people would be in a similar situation), and the balance is not very good when displayed together with sentences and chords. {repr} is a subpackage that is installed with IRkernerl, but you can use this feature to specify options as above to make the figure of any size.
In addition, although this note does not directly execute the code, the BUGS code is also described. (JAGS + {rjags} is a specification that prepares the BUGS code in a separate file and executes the calculation.)
I felt that it is a great advantage to use Jupyter Notebook to be able to combine the R code that summarizes the whole, the BUGS code that defines the simulation model, the calculation result output, and the information that tends to be scattered separately.
The transition of the annual ring width of Sugi in Kyoto calculated by the local level model in this way is illustrated as follows.
Fig. 'NenRin' growth width (by JAGS)
The polygonal line with a marker is the observed value, the thick solid line is the value estimated by MCMC calculation, and the gray area is the 95% confidence interval. (The R code can be obtained from the "Iwanami Data Science" support site.)
Confirmation by {dlm} package
JAGS is an implementation system of MCMC (Marcov chain Monte Carlo). In addition to MCMC, there is an implementation of the Kalman filter system as a calculation method that can handle state space models (dynamic linear models), but in "R", the {dlm} and {KFAS} packages are often used. This time I tried the {dlm} package.
The difference between MCMC and Kalman filter is explained properly in "Iwanami Data Science". (I will quote ...) "The difference between the Kalman filter and MCMC is that the peripheral posterior distribution can be calculated for the state $ x_i $, but it looks like $ s ^ 2 $, $ \ sigma ^ 2 $. For the parameters of the model, the estimation (point estimation) is performed ignoring the spread of the posterior distribution. ”The spread of the parameters of the variance relation, so to speak, the“ variation of variation ”is not very interested (personally). ， This point is not a problem at first.
Below, the R code is quoted from the created Notebook.
## Year and Growth width data(mm)
X <- 1961:1990
Y <- c(4.71, 7.70, 7.97, 8.35, 5.70,
7.33, 3.10, 4.98, 3.75, 3.35,
1.84, 3.28, 2.77, 2.72, 2.54,
3.23, 2.45, 1.90, 2.56, 2.12,
1.78, 3.18, 2.64, 1.86, 1.69,
0.81, 1.02, 1.40, 1.31, 1.57)
tsData <- ts(Y, frequency=1, start=c(1961, 1)) # Time Series data structure
n <- length(tsData)
library(dlm)
# Modeling function and initial parameters
mkModel1 <- function(parm){
dlmModPoly(order = 1, dV=parm[1], dW=parm[2])
}
myparm1 <- c(0.001, 0.001)
# Parameter calculation by fitting
dlmFit1 <- dlmMLE(tsData,
parm = myparm1,
build = mkModel1
)
# Filtering
oFitted1 <- mkModel1(dlmFit1$par)
oFiltered1 <- dlmFilter(tsData, oFitted1)
print(names(oFiltered1))
print(oFiltered1$f)
# Smoothing
oSmoothed1 <- dlmSmooth(oFiltered1)
print(names(oSmoothed1))
print(oSmoothed1$s)
With the above, the calculation result of the Kalman filter can be obtained. Plot this result with {ggplot2}.
library(ggplot2)
base.family <- ""
ribbon.alpha <- 0.3
ribbon.fill <- "black"
alpha.mean <- dropFirst(oSmoothed1$s)
var_p <- unlist(dlmSvd2var(oSmoothed1$U.S, oSmoothed1$D.S))
alpha.up <- dropFirst(oSmoothed1$s) + qnorm(0.025, sd = sqrt(var_p[-1]))
alpha.low <- dropFirst(oSmoothed1$s) + qnorm(0.975, sd = sqrt(var_p[-1]))
p3 <- ggplot(data.frame(Year = X, Width = Y,
Width.mean = alpha.mean,
Width.up = alpha.up,
Width.low = alpha.low)) +
geom_ribbon(aes(x = Year, ymax = Width.up, ymin = Width.low),
fill = ribbon.fill, alpha = ribbon.alpha) +
geom_line(aes(x = Year, y = Width.mean), size = 1) +
geom_line(aes(x = Year, y = Width)) +
geom_point(aes(x = Year, y = Width), size = 2) +
xlab("year") +
ylab("growth width (mm)") +
ylim(0, 10) +
theme_classic(base_family = base.family)
print(p3)
Fig. 'NenRin' growth width (by {dlm})
A closer look reveals a slight difference in the size of the gray area, which indicates the 95% confidence interval, but is almost in agreement with the JAGS results shown above.
Since I usually use Python, I would like to investigate the implementation of Python's Kalman filter system. I tried using Jupyter Notebook this time, and realized that it gives a good view to the data analysis of "R". By the way, "Kalman and Bayesian Filters in Python" is published on Github as an example of Jupyter Notebook, so it is recommended for those who are interested. (Since it is a fairly large amount of content, I am reading it little by little. Animation etc. are also used, and since it is a Jupyter Notebook, you can move the code directly. The link is described below.)
References (web site)
--"Iwanami Data Science" support page https://sites.google.com/site/iwanamidatascience/
- IRkernel - GitHub https://github.com/IRkernel/IRkernel --Local level model --Logic of Blue http://logics-of-blue.com/%E3%83%AD%E3%83%BC%E3%82%AB%E3%83%AB%E3%83%AC%E3%83%99%E3%83%AB%E3%83%A2%E3%83%87%E3%83%AB/
--Reproduce "Introduction to State-Space Time Series Analysis" with R http://elsur.jpn.org/ck/ --{dlm} Manual attached to the package (dlm: an R package for Bayesian analysis of Dynamic Linear Models)
- "Kalman and Bayesian Filters in Python” https://github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python
Recommended Posts