[18846] in s-news-athena
[S] summary: kalman filter in splus
daemon@ATHENA.MIT.EDU (Verena RESSEL)
Wed Aug 4 07:47:17 1999
Message-Id: <199908041130.NAA21354@mail.boku.ac.at>
From: "Verena RESSEL" <vressel@edv2.boku.ac.at>
To: s-news@wubios.wustl.edu
Date: Wed, 4 Aug 1999 13:43:24 +0100
Mime-Version: 1.0
Content-Type: text/enriched; charset=ISO-8859-1
Content-Transfer-Encoding: Quoted-printable
<color><param>0100,0100,0100</param>Thanks to everyone who replied. Here i=
s a summary of the answers
I received.
Verena
-----------------------------------------------------------
<color><param>0000,0000,0000</param>Mrigesh Kshatriya <<<underline><colo=
r><param>0000,8000,0000</param>mkshatriya@zoology.up.ac.za </underline>> <=
color><param>0000,0000,0000</param>wrote:<underline><color><param>0000,800=
0,0000</param>
</underline><color><param>0000,0000,0000</param>I found a time series pack=
age at
<underline><color><param>0000,8000,0000</param>http://anson.ucdavis.edu/~s=
humway/tsa.html</underline><color><param>0000,0000,0000</param> which seem=
s to apply
kalman filter for S.S. analysis.
-----------------------------------------------------------
Don McKenzie <<<underline><color><param>0000,8000,0000</param>dmck@u.washi=
ngton.edu</underline>> <color><param>0000,0000,0000</param>wrote:
There is a posting on statlib that uses EM for the variance estimates.
Also see below for a helpful reply I got to this question (I haven't
done anything about it yet).
---------- Forwarded message ----- Date: Fri, 21 May 1999 09:21:33
From: Paul Gilbert <<<underline><color><param>0000,8000,0000</param>pgilbe=
rt@bank-banque-canada.ca</underline><color><param>0000,0000,0000</param>> =
To: Jorge de la Vega Gongora <<<underline><color><param>0000,8000,0000</pa=
ram>jvega@stat.umn.edu</underline><color><param>0000,0000,0000</param>>
Subject: Re: [S] State space models
There is an old version of my DSE library on Statlib, but you can
find
more recent versions at <<<underline><color><param>0000,8000,0000</param>h=
ttp://www.bank-banque-
canada.ca/pgilbert</underline><color><param>0000,0000,0000</param>> along =
with a User Guide and an HTML version
of the documentation. The library handles multivariate state space
models. The library is not yet working properly with Splus 5. I work
in Unix and am using R and Splus 3.3, though it should work fine
with 3.4. To improve speed the library calls fortran routines for some
calculations, so for reasonable speed you need to be able to
compile fortran, but in theory everything works without the fortran.
I have not implemented EM, but it could be added fairly easily if you
are interested. In general I don't recommend hill climbing
procedures for multivariate state space model estimation.
-------------------------------------------------------------
Paul Gilbert <<<underline><color><param>0000,8000,0000</param>pgilbert@ba=
nk-banque-canada.ca</underline><color><param>0000,0000,0000</param>> wrot=
e:
Below is a copy of the Kalman filter routine clipped from a much
larger time series library I have available at
<<<underline><color><param>0000,8000,0000</param>http://www.bank-banque-ca=
nada.ca/pgilbert</underline><color><param>0000,0000,0000</param>>. One par=
t of this
function sets up parameters to call compiled code which is faster (if
compiled=3D=3DT). You can delete this and use the S code as your
model. It is set up to be called from some other functions, but I think
you will be able to figure out what you need. It also does a few other
things you are probably not interested in doing. I hope it is helpful,
but it will require a lot of changes in order to use for adaptive
estimation.
Please let me know if you get something that works well.
Paul Gilbert
______
KF <<- function (model, data, sampleT =3D NULL, predictT =3D NULL,
error.weights =3D 0, return.state =3D F, return.track =3D F, result =3D N=
ULL,
compiled =3D .DSECOMPILED, warn =3D T, return.debug.info =3D F)
# end of KF
{
# ref. B.D.O.Anderson & J.B.Moore "Optimal Filtering" p.39,44.
# sampleT is the length of data which should be used for
calculating
# one step ahead predictions. y must be at least as
# long as sampleT. If predictT is large than sampleT then the
model is simulated to
# predictT. y is used if it is long enough. u must be at least as
long
as predictT.
# The default result=3D0 returns a list of all the results. Otherwise
only the
# indicated list element is return (eg. result=3D1 return the
likelihood and
# result=3D3 returns the one step ahead predictions.
# Calculate the state, residuals, and likelihood value for the
model:
#
# z(t) =3D Fz(t-1) + Gu(t) + Qe(t)
# y(t) =3D Hz(t) + Rw(t)
#
# or the innovations model:
# z(t) =3D Fz(t-1) + Gu(t)+ Kw(t-1)
# y(t) =3D Hz(t) + w(t)
#
# FF (nxn) is the state transition matrix F.
# H (pxn)is the output matrix H.
# Q (nxn) is the input matrix of the system noise and the noise is
assumed to be white.
# Some authors (eg. Harvey) modify this as rt*qt*rt' where rt is
the matrix for the
# system noise and qt is the noise cov, but that is redundant.
# R (pxp) is the input matrix of the output (measurement) noise,
which is assumed white.
# probably need R if p>n ??
# G (nxp)is the control (input) matrix.
# K (nxp)is the Kalman gain.
# y is the p dimensional output data.
# u is the m dimensional exogenous (input) data.
# z is the n dimensional (estimated) state at time t, E[z(t)|y(t-1), =
u(t)] denoted E[z(t)|t-1].
# Note: In the case where there is no input u this corresponds
to
what
# would usually be called the predicted state - not the
filtered state.
# state is the history of the state.
# Om is the estimated output cov matrix.
# vt is the prediction error.
# pred is the history of the one-step ahead predictions, E[y(t)|y(t-
1),u(t)] denoted E[y(t)|t-1].
# The history of the prediction error is given by y-
pred[1:predictT,]or y-pred[1:sampleT,]
# If error.weights is greater than zero then weighted prediction
# errors are calculated up to the horizon indicated
# by the length of error.weights. The weights are applied to the
squared
# error at each period ahead.
# P is the one step ahead estimate of the state tracking error
matrix at each
# period. Cov{z(t)-E[z(t)|t-1]}
# trackError is the history of P.
# Tracking error pt can only be calculated if Q and R are
provided (gain FALSE).
# Using the Kalman gain K directly these are not necessary
# for the likelihood calculation,
# but the tracking error cannot be calculated.
# If z0 is supplied it is used as the estimate of the state at time 0.=
# If not supplied it is set to zero. # If P0 is supplied it is used as=
the initial tracking error P(t=3D1|t=3D0). # If not supplied it is set to=
I.
# could check that Q is symmetric or positive definite but ...
if (!is.TSm.or.em(model))
TS.error.exit()
if (is.TSestModel(model))
model <<- model$model
if (!is.SS(model))
TS.error.exit(clss =3D "SS")
data <<- TSdata(data)
if (!check.consistent.dimensions(model, data))
stop("dimension error\n")
if (is.null(sampleT))
sampleT <<- output.periods(data)
if (is.null(predictT))
predictT <<- sampleT
if (sampleT > predictT)
stop("sampleT cannot exceed predictT.\n")
if (sampleT > output.periods(data))
stop("sampleT cannot exceed length of data.\n")
if (0 !=3D input.dimension(data)) {
if (input.periods(data) << predictT)
stop("input data must be at least as long as requested
prediction.\n")
if (any(is.na(input.data(data))))
stop("input data cannot contain NAs.\n")
}
if (any(is.na(output.data(data))))
stop("output data cannot contain NAs.\n")
gain <<- is.innov.SS(model)
if (gain & return.track)
warning("Tracking error is zero for an innovations model. track
will not be calculated.")
FF <<- model$F
H <<- model$H
n <<- dim(FF)[2]
p <<- dim(H)[1]
if (is.null(model$G)) {
m <<- 0
# can't call compiled with 0 length arrays
G <<- matrix(0, n, 1)
u <<- matrix(0, predictT, 1)
}
else {
m <<- dim(model$G)[2]
G <<- model$G
u <<- input.data(data)
}
# K or Q,R can be NUll in model, which messes up compiled
if (gain) {
K <<- model$K
#not used
Q <<- matrix(0, 1, 1)
#not used
R <<- matrix(0, 1, 1)
#not used
track <<- array(0, c(1, 1, 1))
}
else {
Q <<- model$Q
# Q assumed square in compiled
if (ncol(Q) << n)
Q <<- cbind(Q, matrix(0, n, n - ncol(Q)))
R <<- model$R
# this is used
K <<- matrix(0, n, p)
#not used
if (return.track)
track <<- array(0, c(predictT, n, n))
else track <<- array(0, c(1, 1, 1))
storage.mode(track) <<- "double"
}
#not used
if (return.state | return.debug.info)
state <<- matrix(0, predictT, n)
else state <<- matrix(0, 1, 1)
storage.mode(state) <<- "double"
# initial state
if (is.null(model$z0))
z <<- rep(0, n)
else z <<- model$z0
# initial state tracking error
# this is not used in innov. models
if (is.null(model$P0))
P <<- diag(1, n)
else P <<- model$P0
# end of S version
if (compiled) {
r <<- .Fortran("kf", pred =3D matrix(0, predictT, p),
as.integer(length(error.weights)),
weighted.sqerror =3D matrix(0, sampleT, p), error.weights =3D
as.double(error.weights),
as.logical(return.state), state =3D state,
as.logical(return.track &
!gain), track =3D track, as.integer(m), as.integer(n),
as.integer(p), sampleT =3D as.integer(sampleT), predictT =3D
as.integer(predictT),
as.integer(output.periods(data)), as.double(u),
as.double(output.data(data)),
as.double(FF), as.double(G), as.double(H), as.double(K),
as.double(Q), as.double(R), as.logical(gain), as.double(z),
as.double(P))[c("pred", "state", "track", "weighted.sqerror")]=
if (all(0 =3D=3D error.weights))
r$weighted.sqerror <<- NULL
}
else {
y <<- output.data(data)
# initial prediction error
vt <<- rep(0, p)
pred <<- matrix(0, predictT, p)
wt.err <<- NULL
if (1 << length(error.weights))
wt.err <<- matrix(0, predictT, p)
if (!gain) {
RR <<- R %*% t(R)
QQ <<- Q %*% t(Q)
}
for (Time in 1:sampleT) {
if (!gain) {
PH <<- P %*% t(H)
ft <<- (H %*% PH) + RR
# force ft to be symmetric
ft <<- (ft + t(ft))/2
K <<- t(solve(ft, t(FF %*% PH)))
# P(t|t-1)
P <<- (FF %*% P %*% t(FF)) - (K %*% H %*% P %*%
t(FF)) + QQ
# force symmetry (eliminate rounding error problems)
P <<- (P + t(P))/2
# P(t|t-1)
# note P(t|t) =3DP -
P%*%t(H)%*%solve(H%*/home/mfa5/gilp/dse/my/src/SCCS/s.dse1b
.hs*%t(H)+RR)%*07/16/99*%P
if (return.track)
track[Time, , ] <<- P
}
# E[z(t)| t-1 ]
z <<- c(FF %*% z) + c(K %*% vt)
if (m !=3D 0)
z <<- z + c(G %*% u[Time, ])
if (return.state | return.debug.info)
state[Time, ] <<- z
# predicted output
pred[Time, ] <<- Ey <<- c(H %*% z)
# prediction error
vt <<- y[Time, ] - Ey
if (any(0 !=3D error.weights)) {
# weighted sq prediction error
wt.err[Time, ] <<- error.weights[1] * vt^2
if (length(error.weights) > 1) {
zh <<- z
for (h in 2:length(error.weights)) if ((Time +
h - 1) <<=3D sampleT) {
zh <<- c(FF %*% zh)
# vt is 0 for h>2
if (h =3D=3D 2)
zh <<- zh + c(K %*% vt)
if (m !=3D 0)
zh <<- zh + c(G %*% u[Time + h - 1, ])
wt.err[Time, ] <<- wt.err[Time, ] + error.weights[h] *=
(y[Time + h - 1, ] - c(H %*% zh))^2
}
}
}
}
# prederror <<- y[1:sampleT,,drop=3DF]-
pred[1:sampleT,,drop=3DF] #
now multi-step prediction to predictT # This requires u but not
y
(y is ignored if it is supplied) if (predictT > sampleT) {
for (Time in (sampleT + 1):predictT) {
z <<- c(FF %*% z)
if (m !=3D 0)
z <<- z + c(G %*% u[Time, ])
if (Time =3D=3D sampleT + 1)
z <<- z + c(K %*% vt)
if (return.state)
state[Time, ] <<- z
# predicted output
pred[Time, ] <<- c(H %*% z)
}
}
r <<- list(pred =3D pred, state =3D state, track =3D track,
weighted.sqerror =3D wt.err)
}
tf <<- truncate.tframe(tframe(output.data(data)), end =3D predictT)
tframe(r$pred) <<- tf
if (!is.null(r$weighted.sqerror))
tframe(r$weighted.sqerror) <<- tf
filter <<- NULL
if (return.state | return.track) {
if (gain | (!return.track))
filter$track <<- NULL
else {
filter$track <<- r$track
tframe(filter$track) <<- tf
}
if (return.state) {
filter$state <<- r$state
tframe(filter$state) <<- tf
}
}
if ((!is.null(result)) && (result =3D=3D "pred"))
return(r$pred)
r <<- append(residual.stats(r$pred, output.data(data), sampleT,
warn =3D warn), list(error.weights =3D error.weights,
weighted.sqerror
=3D
r$weighted.sqerror))
if (return.debug.info)
r$debug.info <<- list(m =3D m, p =3D p, a =3D a, b =3D b, c =3D c,=
G =3D G, FF =3D FF, K =3D K, P =3D P, H =3D H, u =3D u, y =3D =
y,
pred =3D pred, wt.err =3D wt.err, error.weights =3D error.weig=
hts,
sampleT =3D sampleT)
if (is.null(result)) {
r <<- list(estimates =3D r, data =3D data, model =3D model,
filter =3D filter)
return(classed(r, "TSestModel"))
}
else {
if (result =3D=3D "like")
return(r$like[1])
else {
# neg.log.like. from residual.stats
return(r[[result]])
}
}
"should never get to here in KF"
}
---------------------------------------------------------------------
Fernando TUSELL<< <underline><color><param>0000,8000,0000</param>etptupaf@=
bs.ehu.es</underline><color><param>0000,0000,0000</param>> wrote:
You may find useful the following code. I got it from somewhere, a
reference is at the end.
ft.
---------------- code cited in message -------------------
(1) This submission provides a pair of ASCII files to be "sourced"
into S. The resulting S functions, named filter and forecast,
implement Kalman filter recursions using the EM algorithm.
Reference: "An approach to time series smoothing
and forecasting using the EM algorithm."
R. H. Shumway and D. S. Stoffer
J. Time Series Anal. Vol. 3, No. 4 253-264
(2) The information I would like to share:
filter <<-
function(numobs,y,A,Sigma0,mu0,Phi0,Q0,R0,tol=3D0.01,max.iter=3D20,v
erbose=3DT){
# y[,t] rx1
# A[[t]] rxp
# Sigma0 pxp
# mu0 px1
# Phi0 pxp
# Q0 pxp
# R0 rxr
# tol, log likelihood convergence tolerance
# EM estimation of Kalman filter
# Reference: "An approach to time series smoothing
# and forecasting using the EM algorithm."
# R. H. Shumway and D. S. Stoffer
# J. Time Series Anal. Vol. 3, No. 4 253-264
# y(t) =3D A(t)%*%x(t) + v(t), for t=3D1,...,T
# x(t) =3D Phi%*%x(t-1) + w(t), for t=3D1,...,T
# y rx1 observation vector
# A rxp design matrix
# E[v(t)]=3D0, R =3D E[v(t)%*%t(v(t))] rxr observation noise covariance
# x px1 state vector
# Phi pxp transition matrix
# E[w(t)]=3D0, Q =3D E[w(t)%*%t(w(t))] pxp state noise covariance
# initial vector x0, E[x0] =3D mu, Sigma =3D E[(x0-mu)%*%t(x0-mu)]
# Notation:
# sxt =3D E[x(t)|y(1),...,y(s)]
# sPtu =3D E[(x(t)-sxt)%*%t(x(u)-sxu)|y(1),...,y(s)]
N <<- numobs
mu <<- mu0
Phi <<- Phi0
Q <<- Q0
R <<- R0
iter <<- 0
cvg <<- 1
while((iter <<=3D max.iter) && (cvg > tol)){
# gain K pxr and covariances sPtt, tPtt
sPtt <<- vector("list",N)
K <<- vector("list",N)
tPtt <<- vector("list",N+1)
tPtt[[1]] <<- Sigma0
sxt <<- vector("list",N)
txt <<- vector("list",N+1)
txt[[1]] <<- mu
eps <<- vector("list",N)
Sigma <<- vector("list",N)
qf <<- seq(N)*0
det <<- seq(N)*0
# forward recursions
for(n in seq(N)){
sPtt[[n]] <<- Phi%*%tPtt[[n]]%*%t(Phi) + Q
K[[n]] <<-
sPtt[[n]]%*%t(A[[n]])%*%solve(A[[n]]%*%sPtt[[n]]%*%t(A[[n]])+R)
tPtt[[n+1]] <<- sPtt[[n]] - K[[n]]%*%A[[n]]%*%sPtt[[n]]
# filter estimators sxt, txt
sxt[[n]] <<- Phi%*%txt[[n]]
txt[[n+1]] <<- sxt[[n]] + (K[[n]]%*%(y[,n]-(A[[n]]%*%sxt[[n]])))
# innovations eps
eps[[n]] <<- y[,n] - A[[n]]%*%sxt[[n]]
Sgm <<- A[[n]]%*%sPtt[[n]]%*%t(A[[n]]) + R
Sigma[[n]] <<- (t(Sgm)+Sgm)/2
# To ensure its symmetric!
qf[n] <<- t(eps[[n]])%*%solve(Sigma[[n]])%*%eps[[n]]
det[n] <<- prod(eigen( Sigma[[n]] )$values)
# Det for real-valued matrices!
# det[n] <<- prod(diag(chol( Sigma[[n]] )))^2
}
# backward recursions
J <<- vector("list",N)
Txs <<- vector("list",N+1)
Txs[[N+1]] <<- txt[[N+1]]
TPss <<- vector("list",N+1)
TPss[[N+1]] <<- tPtt[[N+1]]
st0 <<- vector("list",N)
ss0 <<- vector("list",N)
St0 <<-
matrix(0,nrow=3Ddim(TPss[[N+1]])[1],ncol=3Ddim(TPss[[N+1]])[2]) Ss0
<<- matrix(0,nrow=3Ddim(TPss[[N+1]])[1],ncol=3Ddim(TPss[[N+1]])[2])
# smoother estimators Txs
for(n in seq(from=3DN,to=3D1)){
J[[n]] <<- tPtt[[n]]%*%t(Phi)%*%solve(sPtt[[n]])
Txs[[n]] <<- txt[[n]] + J[[n]]%*%(Txs[[n+1]]-sxt[[n]])
# covariance TPss
TPss[[n]] <<- tPtt[[n]] +
J[[n]]%*%(TPss[[n+1]]-sPtt[[n]])%*%t(J[[n]])
# Q(,) calculations
st0[[n]] <<- TPss[[n+1]] + Txs[[n+1]]%*%t(Txs[[n+1]])
ss0[[n]] <<- TPss[[n]] + Txs[[n]]%*%t(Txs[[n]])
St0 <<- St0 + st0[[n]]
Ss0 <<- Ss0 + ss0[[n]]
}
loglik <<- -.5*sum(log(det)) - .5*sum(qf)
if(-2*loglik << 0){
out <<- list(mu, Phi, Q, R)
names(out) <<- c("mu", "Phi", "Q", "R")
warning("-2*loglik << 0")
return(out)
}
if(verbose=3D=3DT) print(-2*loglik)
# EM algorithm
# covariance TPsr
TPsr <<- vector("list",N)
p <<- dim(Phi0)[1]
TPsr[[N]] <<- (diag(p) - K[[N]]%*%A[[N]])%*%Phi%*%tPtt[[N]]
for(n in seq(from=3DN,to=3D2)){
TPsr[[n-1]] <<- tPtt[[n]]%*%t(J[[n-1]]) + J[[n]]%*%(TPsr[[n]] -
Phi%*%tPtt[[n]])%*%t(J[[n-1]]) }
# Q(,) calculations
st1 <<- vector("list",N)
St1 <<- matrix(0,nrow=3Ddim(TPsr[[N]])[1],ncol=3Ddim(TPsr[[N]])[2])
rcov <<- matrix(0,nrow=3Ddim(R)[1],ncol=3Ddim(R)[2])
for(n in seq(from=3DN,to=3D1)){
st1[[n]] <<- TPsr[[n]] + Txs[[n+1]]%*%t(Txs[[n]])
St1 <<- St1 + st1[[n]]
rcov <<- rcov +
(y[,n]-A[[n]]%*%Txs[[n+1]])%*%t(y[,n]-A[[n]]%*%Txs[[n+1]]) +
A[[n]]%*%TPss[[n+1]]%*%t(A[[n]]) }
# update mu, Phi, Q, R
# NB: structured parameter matrices can be used by fixing the
appropriate
# elements
mu <<- Txs[[1]]
Phi <<- St1%*%solve(Ss0)
Q <<- matrix(1/N,nrow=3Ddim(Q)[1],ncol=3Ddim(Q)[2])*(St0 -
St1%*%solve(Ss0)%*%t(St1)) R <<-
matrix(1/N,nrow=3Ddim(R)[1],ncol=3Ddim(R)[2])*rcov
iter <<- iter + 1
if(iter>1) cvg <<- (-2*old.loglik) - (-2*loglik)
old.loglik <<- loglik
}
# end while loop
r <<- dim(R0)[1]
BIC <<- loglik - 0.5*(p + p^2 + r*(r+1)/2 + r*(r+1)/2)*log(N) # mu,
Phi, Q, R parameters out <<- list(mu,Phi,Q,R,BIC) names(out) <<-
c("mu","Phi","Q","R","BIC") out
}
----------
forecast <<- function(ahead,numobs,y,A,Sigma0,mu0,Phi0,Q0,R0){
# ahead - the desired number of forecasts of y[,t], t>T
# y[,t] rx1
# A[[t]] rxp
# Sigma0 pxp
# mu0 px1
# Phi0 pxp
# Q0 pxp
# R0 rxr
# EM estimation of Kalman filter
# Reference: "An approach to time series smoothing
# and forecasting using the EM algorithm."
# R. H. Shumway and D. S. Stoffer
# J. Time Series Anal. Vol. 3, No. 4 253-264
# y(t) =3D A(t)%*%x(t) + v(t), for t=3D1,...,T
# x(t) =3D Phi%*%x(t-1) + w(t), for t=3D1,...,T
# y rx1 observation vector
# A rxp design matrix
# E[v(t)]=3D0, R =3D E[v(t)%*%t(v(t))] rxr observation noise covariance
# x px1 state vector
# Phi pxp transition matrix
# E[w(t)]=3D0, Q =3D E[w(t)%*%t(w(t))] pxp state noise covariance
# initial vector x0, E[x0] =3D mu, Sigma =3D E[(x0-mu)%*%t(x0-mu)]
# Notation:
# sxt =3D E[x(t)|y(1),...,y(s)]
# sPtu =3D E[(x(t)-sxt)%*%t(x(u)-sxu)|y(1),...,y(s)]
N <<- numobs + ahead
mu <<- mu0
Phi <<- Phi0
Q <<- Q0
R <<- R0
r <<- dim(A[[1]])[1]
p <<- dim(A[[1]])[2]
# gain K pxr and covariances sPtt, tPtt
sPtt <<- vector("list",N)
K <<- vector("list",N)
tPtt <<- vector("list",N+1)
tPtt[[1]] <<- Sigma0
sxt <<- vector("list",N)
txt <<- vector("list",N+1)
txt[[1]] <<- mu
# forward recursions
for(n in seq(N)){
sPtt[[n]] <<- Phi%*%tPtt[[n]]%*%t(Phi) + Q
# K(t)=3D0, for t>T
if(n<<=3Dnumobs) { K[[n]] <<-
sPtt[[n]]%*%t(A[[n]])%*%solve(A[[n]]%*%sPtt[[n]]%*%t(A[[n]])+R)
}
else { K[[n]] <<- matrix(0,nrow=3Dp,ncol=3Dr) } tPtt[[n+1]] <<-
sPtt[[n]] - K[[n]]%*%A[[n]]%*%sPtt[[n]]
# filter estimators sxt, txt
sxt[[n]] <<- Phi%*%txt[[n]]
if(n<<=3Dnumobs) { txt[[n+1]] <<- sxt[[n]] +
(K[[n]]%*%(y[,n]-(A[[n]]%*%sxt[[n]]))) } else { txt[[n+1]] <<-
sxt[[n]] } }
# backward recursions
J <<- vector("list",N)
Txs <<- vector("list",N+1)
Txs[[N+1]] <<- txt[[N+1]]
TPss <<- vector("list",N+1)
TPss[[N+1]] <<- tPtt[[N+1]]
# smoother estimators Txs
for(n in seq(from=3DN,to=3D1)){
J[[n]] <<- tPtt[[n]]%*%t(Phi)%*%solve(sPtt[[n]])
Txs[[n]] <<- txt[[n]] + J[[n]]%*%(Txs[[n+1]]-sxt[[n]])
# covariance TPss
TPss[[n]] <<- tPtt[[n]] +
J[[n]]%*%(TPss[[n+1]]-sPtt[[n]])%*%t(J[[n]]) }
# covariance TPsr
TPsr <<- vector("list",N)
TPsr[[N]] <<- (diag(p) - K[[N]]%*%A[[N]])%*%Phi%*%tPtt[[N]]
for(n in seq(from=3DN,to=3D2)){
TPsr[[n-1]] <<- tPtt[[n]]%*%t(J[[n-1]]) + J[[n]]%*%(TPsr[[n]] -
Phi%*%tPtt[[n]])%*%t(J[[n-1]]) }
Ty <<- matrix(0,nrow=3Dr,ncol=3DN)
Tx <<- matrix(0,nrow=3Dp,ncol=3DN)
for(i in seq(from=3DN+1,to=3D2)){
Ty[,i-1] <<- A[[i-1]]%*%Txs[[i]]
Tx[,i-1] <<- Txs[[i]]
}
out <<- list(Ty,Tx)
names(out) <<- c("obs.pred","state.pred")
out
}
----------
(3) GENERAL DISCLAIMER
This program is free software; you can redistribute it
and/or modify it under the terms of the GNU General Public
License as published by the Free Software Foundation; either
version 1, or (at your option) any later version.
This program is distributed in the hope that it will be
useful, but WITHOUT ANY WARRANTY; without even the
implied
warranty of MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE. See the GNU General Public License
for more details.
For a copy of the GNU General Public License write to the
Free Software Foundation, Inc., 675 Mass Ave, Cambridge,
MA 02139, USA. In short: you may use this code any way you
like, as long as you don't charge money for it, remove this notice, or
hold anyone liable for its results.
Gary Nan Tie
(2i2) 902-7265
<underline><color><param>0000,8000,0000</param>nantig@fi.gs.com
</underline>----------------------------------------------<underline>
<nofill>
------------------------------------------------------------------------
DI Verena Ressel ph: +43 1 36006 5525
email: vressel@edv2.boku.ac.at
Institut f=FCr Wasserwirtschaft, Hydrologie
und konstruktiven Wasserbau (IWHW)
Universit=E4t f=FCr Bodenkultur Wien
Muthgasse 18
A-1190 Wien
AUSTRIA
------------------------------------------------------------------------
-----------------------------------------------------------------------
This message was distributed by s-news@wubios.wustl.edu. To unsubscribe
send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
message: unsubscribe s-news
