没有合适的资源?快使用搜索试试~ 我知道了~
首页kalman filtering in R(R语言实现卡尔曼滤波)
kalman filtering in R(R语言实现卡尔曼滤波)
4星 · 超过85%的资源 需积分: 49 102 下载量 131 浏览量
更新于2023-03-16
评论 4
收藏 842KB PDF 举报
卡尔曼滤波(平方根算法、序贯算法、线性卡尔曼滤波、扩展卡尔曼滤波、区间卡尔曼滤波、系统误差和观测误差有关的卡尔曼滤波)在R语言中的实现,详细讲述了原理和其中工具包的使用,英文原版。
资源详情
资源评论
资源推荐
JSS
Journal of Statistical Software
March 2011, Volume 39, Issue 2. http://www.jstatsoft.org/
Kalman Filtering in R
Fernando Tusell
University of the Basque Country
Abstract
Support in R for state space estimation via Kalman filtering was limited to one package,
until fairly recently. In the last five years, the situation has changed with no less than
four additional packages offering general implementations of the Kalman filter, including in
some cases smoothing, simulation smoothing and other functionality. This paper reviews
some of the offerings in R to help the prospective user to make an informed choice.
Keywords: state space models, Kalman filter, time series, R.
1. Introduction
The Kalman filter is an important algorithm, for which relatively little support existed in R
(R Development Core Team 2010) up until fairly recently. Perhaps one of the reasons is the
(deceptive) simplicity of the algorithm, which makes it easy for any prospective user to throw
in his/her own quick implementation.
While direct transcription of the equations of the Kalman filter as they appear in many
engineering or time series books may be sufficient for some applications, an all-around im-
plementation requires more complex coding. In Section 2, we provide a short overview of
available algorithms. It is against this background that we discuss in Section 3 the particular
choices made by four packages offering fairly general Kalman filtering in R, with some mention
of functions in other packages which cater to particular needs.
Kalman filtering is a large topic. In the sequel we focus on linear Gaussian models and
their estimation, which is what the packages we review offer in common (and the foundation
on which most anything else rests). Other functionalities present in some of the packages
examined include filtering and estimation of non-Gaussian models, simulation and disturbance
smoothing and functions to help with the Bayesian analysis of dynamic linear models, etc.
none of which are assessed.
2 Kalman Filtering in R
2. Kalman filter algorithms
We shall consider a fairly general state-space model specification, sufficient for the purpose
of the discussion to follow in Section 3, even if not the most comprehensive. The notation
follows Harvey (1989). Let
α
t
= c
t
+ T
t
α
t−1
+ R
t
η
t
(1)
y
t
= d
t
+ Z
t
α
t
+
t
(2)
where η
t
∼ N(0, Q
t
) and
t
∼ N(0, H
t
). The state equation (1) describes the dynamics of
the state vector α
t
, driven by deterministic (c
t
) and stochastic (η
t
) inputs. The observation
(or measurement) equation links the observed response y
t
with the unobserved state vector,
with noise
t
and (possibly) deterministic inputs d
t
. Matrices T
t
, R
t
, Z
t
, Q
t
and H
t
may
depend on a vector of parameters, θ, and be time varying or constant over time. The noises η
t
and
t
are assumed serially and also mutually uncorrelated, i.e., E[η
t
>
s
] = 0 for all t, s. The
last assumption and the gaussianity of η
t
and
t
can be dispensed with; see e.g., Anderson
and Moore (1979).
2.1. The Kalman filter.
The Kalman filter equations, with slight notational variations, are standard in any textbook:
see, e.g., Anderson and Moore (1979), Simon (2006), Durbin and Koopman (2001), Grewal
and Andrews (2001), West and Harrison (1997) or Shumway and Stoffer (2006), to name only
a few. We reproduce those equations here, however, as repeated reference is made to them in
the sequel. Define
a
t−1
= E[α
t−1
|y
0
, . . . , y
t−1
] (3)
P
t−1
= E[(α
t−1
− a
t−1
)(α
t−1
− a
t−1
)
>
] ; (4)
estimates of the state vector and its covariance matrix at time t with information available
at time t − 1, a
t|t−1
and P
t|t−1
respectively, are given by the time update equations
a
t|t−1
= T
t
a
t−1
+ c
t
(5)
P
t|t−1
= T
t
P
t−1
T
>
t
+ R
t
Q
t
R
>
t
. (6)
Let F
t
= Z
t
P
t|t−1
Z
t
>
+ H
t
. If a new observation is available at time t, then a
t|t−1
and
P
t|t−1
can be updated with the measurement update equations
a
t
= a
t|t−1
+ P
t|t−1
Z
>
t
F
−1
t
(y
t
− Z
t
a
t|t−1
− d
t
) (7)
P
t
= P
t|t−1
− P
t|t−1
Z
t
F
−1
t
Z
>
t
P
t|t−1
. (8)
Equations (5)–(6) and (7)–(8) taken together make up the Kalman filter. Substituting (5)
in (7) and (6) in (8), a single set of equations linking a
t−1
and P
t−1
to a
t
and P
t
can be
obtained. In order to start the iteration we need initial values of a
−1
and P
−1
(or a
0|−1
and
P
0|−1
).
The filter equations (6) and (8) propagate the covariance matrix of the state, and are said
to define a covariance filter (CF). Equivalent equations can be written which propagate the
matrix P
−1
t
, giving an information filter (IF). Information filters require in general more
Journal of Statistical Software 3
computational effort. One possible advantage is that they provide a natural way to specify
complete uncertainty about the initial value of a component of the state: we can set the
corresponding diagonal term in the information matrix to zero. With a covariance filter, we
have to set the corresponding variance in P
0|−1
to a “large” number, or else use exact diffuse
initialization, an option described below.
Direct transcription of the equations making up the Kalman filter into computer code is
straightforward. It was soon noticed, though, that the resulting programs suffered from
numerical instability; see for instance Bucy and Joseph (1968). In particular, buildup of
floating point errors in equation (8) may eventually yield non symmetric or non positive
definite P
t
matrices. An alternative to equation (8) (more expensive from the computational
point of view) is:
P
t
= (I − K
t
Z
t
)P
t|t−1
(I − K
t
Z
t
)
>
+ K
t
H
t
K
>
t
(9)
with K
t
= P
t|t−1
Z
>
t
F
−1
t
; but even this (“Joseph stabilized form”, see Bucy and Joseph
(1968), p. 175) may prove insufficient to prevent roundoff error degeneracy in the filter. A
detailed reference describing the pitfalls associated with the numerical implementation of the
Kalman filter is Bierman (1977); see also Grewal and Andrews (2001), Chap. 6 and Anderson
and Moore (1979), Chap. 6. Square root algorithms, that we discuss next, go a long way in
improving the numerical stability of the filter.
2.2. Square root algorithms
Consider a matrix S
t
such that P
t
= S
t
S
t
>
; square root covariance filter (SRCF) algorithms
propagate S
t
instead of P
t
with two main benefits (cf. Anderson and Moore (1979), § 6.5):
(i) Re-constitution of P
t
from S
t
will always yield a symmetric non negative matrix, and (ii)
The numerical condition of S
t
will in general be much better than that of P
t
. If instead of
the covariance matrix we choose to factor the information matrix P
−1
t
we have a square root
information filter algorithm (SRIF).
It is easy to produce a replacement for the time update equation (6). Consider an orthogonal
matrix G such that:
M
0
= G
S
>
t−1
T
>
t
(Q
1/2
)
T
R
>
t
(10)
where M is an upper triangular matrix. Matrix G can be constructed in a variety of ways,
including repeated application of Householder or Givens transforms. Multiplication of the
transpose of (10) by itself produces, taking into account the orthogonality of G,
M
>
M = T
t
S
t−1
S
>
t−1
T
>
t
+ R
t
Q
1/2
t
(Q
1/2
t
)
T
R
t
>
(11)
= T
t
P
t−1
T
>
t
+ R
t
Q
t
R
>
t
; (12)
comparison with (6) shows that M is a possible choice for S
t|t−1
. Likewise, it can be shown
(Simon 2006, Section 6.3.4) that the orthogonal matrix G
∗
such that
(H
t
+ Z
t
P
t|t−1
Z
>
t
) K
∗
t
>
0 M
∗
= G
∗
"
(H
1/2
t
)
T
0
S
>
t|t−1
Z
>
t
S
>
t|t−1
#
(13)
produces in the left-hand side of (13) a block M
∗
that can be taken as S
t
, and thus performs
the measurement update.
4 Kalman Filtering in R
Although the matrices performing the block triangularizations described in equations (10)
and (13) can be obtained quite efficiently (Lawson and Hanson 1974; Gentle 2007), clearly
the computational effort is greater than that required by the time and measurement updates
in equations (6) and (8); square root filtering does have a cost.
In the equations above G and G
∗
can be chosen so that M and M
∗
are Cholesky factors
of the corresponding covariance matrices. This needs not be so, and other factorizations
are possible. In particular, the factors in the singular value decomposition of P
t−1
can be
propagated: see for instance Zhang and Li (1996) and Appendix B of Petris et al. (2009). Also,
we may note that time and measurement updates can be merged in a single triangularization
(see Anderson and Moore 1979, p. 148, and Vanbegin and Verhaegen 1989).
2.3. Sequential processing
In the particular case where H
t
is diagonal, the components of y
t
>
= (y
1t
, . . . , y
pt
) are
uncorrelated. We may pretend that we observe one y
it
at a time and perform p univariate
measurement updates similar to (7)–(8), followed by a time update (5)–(6) – see Durbin and
Koopman (2001, Section 6.4) or Anderson and Moore (1979) for details.
The advantage of sequential processing is that F
t
becomes 1 ×1 and the inversion of a p ×p
matrix in equation (7)–(8) is avoided. Clearly, the situation where we stand to gain most
from this strategy is when the dimension of the observation vector y
t
is large.
Sequential processing can be combined with square root covariance and information filters,
although in the latter case the computational advantages seem unclear (Anderson and Moore
1979, p. 142; see also Bierman 1977).
Aside from the case where H
t
is diagonal, sequential processing may also be used when H
t
is block diagonal – in which case we can perform a sequence of reduced dimension updates
– or when it can be reduced to diagonal form by a linear transformation. In the last case,
assuming full rank, let H
−1/2
t
be a square root of H
−1
t
. Multiplying (2) through by H
−1/2
t
we get
y
∗
t
= d
∗
t
+ Z
∗
t
α
t
+
∗
t
(14)
with E[
∗
t
∗
t
>
] = I. If matrix H
t
is time-invariant, the same linear transformation will
decorrelate the measurements at all times.
2.4. Smoothing and the simulation smoother
The algorithms presented produce predicted a
t|t−1
or filtered a
t
values of the state vector
α
t
. Sometimes it is of interest to estimate a
t|N
for 0 < t ≤ N, i.e., E[α
t
|y
0
, . . . , y
N
], the
value of the state vector given all past and future observations. It turns out that this can
be done running once the Kalman filter and then a recursion backwards in time (Durbin and
Koopman 2001, Section 4.3, Harvey 1989, Section 3.6).
In some cases, and notably for the Bayesian analysis of the state space model, it is of interest
to generate random samples of state and disturbance vectors, conditional on the observations
y
0
, . . . , y
N
. Fr
¨
uhwirth-Schnatter (1994) and Carter and Kohn (1994) provided algorithms to
that purpose, improved by de Jong (1995). Durbin and Koopman (2001, Section 4.7) draws
on work of the last author; the algorithm they present is fairly involved. Recently, Durbin
and Koopman (2002) have provided a much simpler algorithm for simulation smoothing of
the Gaussian state space model; see also Strickland et al. (2009).
Journal of Statistical Software 5
2.5. Exact diffuse initial conditions
As mentioned above, the Kalman filter iteration needs starting values a
−1
and P
−1
. When
nothing is known about the initial state, a customary practice has been to set P
−1
with
“large” elements along the main diagonal, to reflect our uncertainty. This practice has been
criticized on the ground that
“While [it can] be useful for approximate exploratory work, it is not recommended
for general use, since it can lead to large rounding errors.”
(Durbin and Koopman 2001, p. 101). Grewal and Andrews (2001, Section 6.3.2) show how
this may come about when elements of P
−1
are set to values “too large” relative to the mea-
surement variances. An alternative is to use an information filter with the initial information
matrix (or some diagonal elements of it) set to zero. Durbin and Koopman (2001, Section 5.3)
advocate a different approach: see also Koopman and Durbin (2003) and Koopman (1997).
The last reference discusses several alternatives and their respective merits.
2.6. Maximum likelihood estimation
An important special case of the state space model is that in which some or all of the matrices
T
t
, R
t
, Z
t
, Q
t
and H
t
in equations (1)–(2) are time invariant and depend on a vector of
parameters, θ that we seek to estimate.
Assuming α
0
∼ N (a
0
, P
0
) with both a
0
and P
0
known, the log likelihood is given by,
L(θ) = log p(y
0
, . . . , y
N
|θ)
= −
(N + 1)p
2
log(2π) −
1
2
N
X
t=0
log |F
t
| + e
>
t
F
−1
t
e
t
(15)
where e
t
= y
t
− Z
t
a
t
. (Durbin and Koopman 2001, Section 7.2; the result goes back to
Schweppe 1965). Except for |F
t
|, all other quantities in (15) are computed when running
the Kalman filter (cf. equations (5)–(8)); therefore, the likelihood is easily computed. (If
we resort to sequential processing, Section 2.3, the analog of equation (15) does not require
computation of determinants nor matrix inversions.)
Maximum likelihood estimation can therefore be implemented easily: it suffices to write a
routine computing (15) as a function of the parameters θ and use R functions such as optim
or nlminb to perform the optimization. An alternative is to use the EM algorithm which quite
naturally adapts to this likelihood maximization problem, but that approach has generally
been found slower than the use of quasi-Newton methods, see for instance Shumway and
Stoffer (2006), p. 345.
3. Kalman filtering in R
There are several packages available from the Comprehensive R Archive Network (CRAN)
offering general Kalman filter capabilities, plus a number of functions scattered in other
packages which cater to special models or problems. We describe five of those packages in
chronological order of first appearance on CRAN.
剩余26页未读,继续阅读
feiyihexin
- 粉丝: 26
- 资源: 1
上传资源 快速赚钱
- 我的内容管理 收起
- 我的资源 快来上传第一个资源
- 我的收益 登录查看自己的收益
- 我的积分 登录查看自己的积分
- 我的C币 登录后查看C币余额
- 我的收藏
- 我的下载
- 下载帮助
会员权益专享
最新资源
- RTL8188FU-Linux-v5.7.4.2-36687.20200602.tar(20765).gz
- c++校园超市商品信息管理系统课程设计说明书(含源代码) (2).pdf
- 建筑供配电系统相关课件.pptx
- 企业管理规章制度及管理模式.doc
- vb打开摄像头.doc
- 云计算-可信计算中认证协议改进方案.pdf
- [详细完整版]单片机编程4.ppt
- c语言常用算法.pdf
- c++经典程序代码大全.pdf
- 单片机数字时钟资料.doc
- 11项目管理前沿1.0.pptx
- 基于ssm的“魅力”繁峙宣传网站的设计与实现论文.doc
- 智慧交通综合解决方案.pptx
- 建筑防潮设计-PowerPointPresentati.pptx
- SPC统计过程控制程序.pptx
- SPC统计方法基础知识.pptx
资源上传下载、课程学习等过程中有任何疑问或建议,欢迎提出宝贵意见哦~我们会及时处理!
点击此处反馈
安全验证
文档复制为VIP权益,开通VIP直接复制
信息提交成功
评论2