[Robkalman-commits] r76 - branches/robKalman_2012/pkg/robKalman/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Jan 16 16:59:10 CET 2015
Author: ruckdeschel
Date: 2015-01-16 16:59:10 +0100 (Fri, 16 Jan 2015)
New Revision: 76
Added:
branches/robKalman_2012/pkg/robKalman/R/recSmooth4.R
Modified:
branches/robKalman_2012/pkg/robKalman/R/allClass.R
branches/robKalman_2012/pkg/robKalman/R/recFilter4.R
branches/robKalman_2012/pkg/robKalman/R/updateinitSSPredOrFilterRet.R
Log:
after session 2015 01 16
Modified: branches/robKalman_2012/pkg/robKalman/R/allClass.R
===================================================================
--- branches/robKalman_2012/pkg/robKalman/R/allClass.R 2014-04-16 16:34:59 UTC (rev 75)
+++ branches/robKalman_2012/pkg/robKalman/R/allClass.R 2015-01-16 15:59:10 UTC (rev 76)
@@ -26,6 +26,9 @@
setClassUnion("OptionalMatrix",
c("NULL", "matrix")
)
+setClassUnion("OptionalArray",
+ c("NULL", "array")
+ )
setClassUnion("OptionalList",
c("list","NULL")
)
@@ -235,6 +238,7 @@
representation = representation(values = "matrix",
call = "OptionalListOfCalls",
variances = "array",
+ lag1variances = "OptionalArray",
uExo = "OptionalMatrix",
wExo = "OptionalMatrix",
dots.propagated = "OptionalList",
@@ -282,16 +286,11 @@
steps = "SSFilterOrSmoother")
)
setClass("SSOutput",
- representation = representation(init.cl = "SSInitialized",
- prep.cl = "OptionalSSPreparedRet",
- pred.cl = "SSPredictedRet",
- filt.cl = "SSFilteredRet",
- smooth.cl = "OptionalSSSmoothedRet",
- init.rob = "OptionalSSInitialized",
- pred.rob = "OptionalSSPredictedRet",
- filt.rob = "OptionalSSFilteredRet",
- smooth.rob = "OptionalSSSmoothedRet",
- prep.rob = "OptionalSSPreparedRet")
+ representation = representation(init = "list",
+ prep = "OptionalList",
+ pred = "list",
+ filt = "list",
+ smooth = "OptionalList")
)
setClass("SSrecResult",
Modified: branches/robKalman_2012/pkg/robKalman/R/recFilter4.R
===================================================================
--- branches/robKalman_2012/pkg/robKalman/R/recFilter4.R 2014-04-16 16:34:59 UTC (rev 75)
+++ branches/robKalman_2012/pkg/robKalman/R/recFilter4.R 2015-01-16 15:59:10 UTC (rev 76)
@@ -1,142 +1,143 @@
-#######################################################
-##
-## recursive filter algorithm for Kalman filter routines, S4
-## author: Bernhard Spangl
-## version: 0.2 (changed: 2013-07-16, created: 2013-05-08)
-##
-#######################################################
-
-# _<
-## ist erledigt durch initSSPredOrFiltRet() in updateinitSSPredOrFiltRet.R
-# _>
-# initPsRet <- function(SSM,tt, exosDim ){
-# ###exosDim ist die Dimensionierung der Rueckgabewerte der exoFct
-#
-# if modell has uexo uexos = matrix(...) else uexos = NULL
-# psret0 <- new("SSPredictedRet",
-# values=matrix(0,Model at pdim,length(tt)+1),
-# call = vector("list",length(tt)+1 ),
-# variance = array(0,dim=c(Model at pdim,Model at pdim,length(tt)+1)),
- # uexo = uexos,...)
-# return(psret0)
-# }
-
-
-recFilter <- function (Model,
- Obs,
- times,
- Steps,
- ...)
-{
- ## Model ... object of S4 class 'SSM'
- ## Obs ... object of S4 class 'SSObs'
- ## times ... object of S4 class 'SStimes'
- ## Steps ... object of S4 class 'SSFilterOrSmoother'
- call <- match.call()
- dots.propagated <- list(...)
-
- ## unwrapping:
- initEq <- Model at initEq
- stateEq <- Model at stateEq
- obsEq <- Model at obsEq
-
- ## time management:
- tt <- times at times
- inX <- times at inX
- tT <- length(tt)+1
- tY <- sum(inX)
- loopIndex <- 1:length(tt)
-
- nrSteps <- length(Steps)
-
- ## initialization of resulting objects:
- ps <- vector("list", nrSteps)
- cs <- vector("list", nrSteps)
- psRet <- vector("list", nrSteps)
- csRet <- vector("list", nrSteps)
-
- withuExo <- !is.null(stateEq at uExofct)
- withwExo <- !is.null(obsEq at wExofct)
- withdots.prop <- (length(dots.propagated)>0)
-
- ### noch herauszufinden: woher findet man raus,
- ## ob der prepstep/predstep/corrstep control/Diagnostic hat..
-
- withcontrol.prep <- ##!is.null(stateEq at uExofct)
- withDiagnostic.prep <- ##!!is.null(stateEq at uExofct)
-
- withcontrol.pred <- ##!!is.null(stateEq at uExofct)
- withDiagnostic.pred <- ##!!is.null(stateEq at uExofct)
-
- withcontrol.corr <- ##!!is.null(stateEq at uExofct)
- withDiagnostic.corr <- ##!!is.null(stateEq at uExofct)
-
- if (prep) prepret0 <- initSSPredOrFiltRet(Model at pdim, tT,
- Model at pdim, tT, Model at qdim, tY,
- withuExo, withwExo, withdots.prop,
- withcontrol.prep, withDiagnostic.prep)
- psret0 <- <- initSSPredOrFiltRet(Model at pdim, tT,
- Model at pdim, tT, Model at qdim, tY,
- withuExo, withwExo, withdots.prop,
- withcontrol.pred, withDiagnostic.pred)
- csret0 <- <- initSSPredOrFiltRet(Model at pdim, tY,
- Model at pdim, tT, Model at qdim, tY,
- withuExo, withwExo, withdots.prop,
- withcontrol.corr, withDiagnostic.corr)
-
- for(iStep in 1:nrStep){
- if (prep) prepRet[[iStep]] <- prepret0
- psRet[[iStep]] <- psret0
- csRet[[iStep]] <- csret0
- }
- ## initialization: iStep = index within different procedures
- for (iStep in 1:nrStep) {
- cs[[iStep]] <- Steps[[iStep]]@initStep(initEq=initEq, ...)
- }
- iniRet <- cs
- ### ab hier weiter wie oben!, 2013-07-16 ###
- ### Schleife über Zeitpunkte t mit entsprechendem Time-Management ###
-
- iy <- 0
- for(ix in loopIndex){
- ## preparation: TBD!
- for (iStep in 1:nrStep) {
- if (!is.null(Steps[[iStep]]@prepStep) {
- ### to be filled
- ##preps[[iStep]] <- Steps[[iStep]]@prepStep(i=ix, t=tt[ix],
- # PredOrFilt=cs[[iStep]],
- # stateEq=stateEq, ...)
- #prepsRet[[iStep]] <- updateSSPredOrFilt(prepsRet[[iStep]], preps[[iStep]],
- # ix)
- #cs[[iStep]] <- preps[[iStep]]
- }
- }
-
- ## prediction:
- for (iStep in 1:nrStep) {
- ps[[iStep]] <- Steps[[iStep]]@predStep(i=ix, t=tt[ix],
- PredOrFilt=cs[[iStep]],
- stateEq=stateEq, ...)
- psRet[[iStep]] <- updateSSPredOrFilt(psRet[[iStep]], ps[[iStep]],
- ix)
- }
-
- ## correction:
- if(inX[ix]){ ## have an observation available
- iy <- iy + 1
- for (iStep in 1:nrStep) {
- cs[[iStep]] <- Steps[[iStep]]@corrStep(i=iy, t=tt[ix], Obs=Obs,
- PredOrFilt=ps[[iStep]],
- obsEq=obsEq, ...)
- csRet[[iStep]] <- updateSSPredOrFilt(csRet[[iStep]],
- cs[[iStep]], iy)
- }
- }else{
- cs <- ps
- }
- }
- zusammenhängen von list(initRet, prepRet, psRet, csRet) ## SSOutput ist jetzt
- eine Liste mit 4 Elementen; jedes einzelne der 4 ist wieder Liste
- und enthält zB ideal robust -> unterschied zu Output.dia
-
+#######################################################
+##
+## recursive filter algorithm for Kalman filter routines, S4
+## author: Bernhard Spangl
+## version: 0.2 (changed: 2013-07-16, created: 2013-05-08)
+##
+#######################################################
+
+# _<
+## ist erledigt durch initSSPredOrFiltRet() in updateinitSSPredOrFiltRet.R
+# _>
+# initPsRet <- function(SSM,tt, exosDim ){
+# ###exosDim ist die Dimensionierung der Rueckgabewerte der exoFct
+#
+# if modell has uexo uexos = matrix(...) else uexos = NULL
+# psret0 <- new("SSPredictedRet",
+# values=matrix(0,Model at pdim,length(tt)+1),
+# call = vector("list",length(tt)+1 ),
+# variance = array(0,dim=c(Model at pdim,Model at pdim,length(tt)+1)),
+ # uexo = uexos,...)
+# return(psret0)
+# }
+
+
+recFilter <- function (Model,
+ Obs,
+ times,
+ Steps,
+ ...)
+{
+ ## Model ... object of S4 class 'SSM'
+ ## Obs ... object of S4 class 'SSObs'
+ ## times ... object of S4 class 'SStimes'
+ ## Steps ... object of S4 class 'SSFilterOrSmoother'
+ call <- match.call()
+ dots.propagated <- list(...)
+
+ ## unwrapping:
+ initEq <- Model at initEq
+ stateEq <- Model at stateEq
+ obsEq <- Model at obsEq
+
+ ## time management:
+ tt <- times at times
+ inX <- times at inX
+ tT <- length(tt)+1
+ tY <- sum(inX)
+ loopIndex <- 1:length(tt)
+
+ nrSteps <- length(Steps)
+
+ ## initialization of resulting objects:
+ ps <- vector("list", nrSteps)
+ cs <- vector("list", nrSteps)
+ psRet <- vector("list", nrSteps)
+ csRet <- vector("list", nrSteps)
+
+ withuExo <- !is.null(stateEq at uExofct)
+ withwExo <- !is.null(obsEq at wExofct)
+ withdots.prop <- (length(dots.propagated)>0)
+
+ ### noch herauszufinden: woher findet man raus,
+ ## ob der prepstep/predstep/corrstep control/Diagnostic hat..
+
+ withcontrol.prep <- ##!is.null(stateEq at uExofct)
+ withDiagnostic.prep <- ##!!is.null(stateEq at uExofct)
+
+ withcontrol.pred <- ##!!is.null(stateEq at uExofct)
+ withDiagnostic.pred <- ##!!is.null(stateEq at uExofct)
+
+ withcontrol.corr <- ##!!is.null(stateEq at uExofct)
+ withDiagnostic.corr <- ##!!is.null(stateEq at uExofct)
+
+ if (prep) prepret0 <- initSSPredOrFiltRet(Model at pdim, tT,
+ Model at pdim, tT, Model at qdim, tY,
+ withuExo, withwExo, withdots.prop,
+ withcontrol.prep, withDiagnostic.prep)
+ psret0 <- initSSPredOrFiltRet(Model at pdim, tT,
+ Model at pdim, tT, Model at qdim, tY,
+ withuExo, withwExo, withdots.prop,
+ withcontrol.pred, withDiagnostic.pred)
+ csret0 <- initSSPredOrFiltRet(Model at pdim, tY,
+ Model at pdim, tT, Model at qdim, tY,
+ withuExo, withwExo, withdots.prop,
+ withcontrol.corr, withDiagnostic.corr)
+
+ for(iStep in 1:nrStep){
+ if (prep) prepRet[[iStep]] <- prepret0
+ psRet[[iStep]] <- psret0
+ csRet[[iStep]] <- csret0
+ }
+ ## initialization: iStep = index within different procedures
+ for (iStep in 1:nrStep) {
+ cs[[iStep]] <- Steps[[iStep]]@initStep(initEq=initEq, ...)
+ }
+ iniRet <- cs
+ ### ab hier weiter wie oben!, 2013-07-16 ###
+ ### Schleife über Zeitpunkte t mit entsprechendem Time-Management ###
+
+ iy <- 0
+ for(ix in loopIndex){
+ ## preparation: TBD!
+ for (iStep in 1:nrStep) {
+ if (!is.null(Steps[[iStep]]@prepStep) {
+ ### to be filled
+ ##preps[[iStep]] <- Steps[[iStep]]@prepStep(i=ix, t=tt[ix],
+ # PredOrFilt=cs[[iStep]],
+ # stateEq=stateEq, ...)
+ #prepsRet[[iStep]] <- updateSSPredOrFilt(prepsRet[[iStep]], preps[[iStep]],
+ # ix)
+ #cs[[iStep]] <- preps[[iStep]]
+ }
+ }
+
+ ## prediction:
+ for (iStep in 1:nrStep) {
+ ps[[iStep]] <- Steps[[iStep]]@predStep(i=ix, t=tt[ix],
+ PredOrFilt=cs[[iStep]],
+ stateEq=stateEq, ...)
+ psRet[[iStep]] <- updateSSPredOrFilt(psRet[[iStep]], ps[[iStep]],
+ ix)
+ }
+
+ ## correction:
+ if(inX[ix]){ ## have an observation available
+ iy <- iy + 1
+ for (iStep in 1:nrStep) {
+ cs[[iStep]] <- Steps[[iStep]]@corrStep(i=iy, t=tt[ix], Obs=Obs,
+ PredOrFilt=ps[[iStep]],
+ obsEq=obsEq, ...)
+ csRet[[iStep]] <- updateSSPredOrFilt(csRet[[iStep]],
+ cs[[iStep]], iy)
+ }
+ }else{
+ cs <- ps
+ }
+ }
+### to do
+# zusammenhaengen von list(initRet, prepRet, psRet, csRet) ## SSOutput ist jetzt
+# eine Liste mit 4 Elementen; jedes einzelne der 4 ist wieder Liste
+# und enthaelt zB ideal robust -> unterschied zu Output.dia
+
}
\ No newline at end of file
Added: branches/robKalman_2012/pkg/robKalman/R/recSmooth4.R
===================================================================
--- branches/robKalman_2012/pkg/robKalman/R/recSmooth4.R (rev 0)
+++ branches/robKalman_2012/pkg/robKalman/R/recSmooth4.R 2015-01-16 15:59:10 UTC (rev 76)
@@ -0,0 +1,98 @@
+#######################################################
+##
+## recursive filter algorithm for Kalman filter routines, S4
+## author: Bernhard Spangl
+## version: 0.2 (changed: 2013-07-16, created: 2013-05-08)
+##
+#######################################################
+
+# _<
+## ist erledigt durch initSSPredOrFiltRet() in updateinitSSPredOrFiltRet.R
+# _>
+# initPsRet <- function(SSM,tt, exosDim ){
+# ###exosDim ist die Dimensionierung der Rueckgabewerte der exoFct
+#
+# if modell has uexo uexos = matrix(...) else uexos = NULL
+# psret0 <- new("SSPredictedRet",
+# values=matrix(0,Model at pdim,length(tt)+1),
+# call = vector("list",length(tt)+1 ),
+# variance = array(0,dim=c(Model at pdim,Model at pdim,length(tt)+1)),
+ # uexo = uexos,...)
+# return(psret0)
+# }
+
+
+recSmoother <- function (Model,
+ Obs,
+ times,
+ Steps,
+ FilterOutput,
+ ...)
+{
+ ## Model ... object of S4 class 'SSM'
+ ## Obs ... object of S4 class 'SSObs'
+ ## times ... object of S4 class 'SStimes'
+ ## Steps ... object of S4 class 'SSFilterOrSmoother'
+ ## FilterOutput ... object of S4 class 'SSOutput'
+ call <- match.call()
+ dots.propagated <- list(...)
+
+ ## unwrapping:
+ initEq <- Model at initEq
+ stateEq <- Model at stateEq
+ obsEq <- Model at obsEq
+
+ ## time management:
+ tt <- times at times
+ inX <- times at inX
+ tT <- length(tt)+1
+ tY <- sum(inX)
+ loopIndex <- 1:length(tt)
+
+ nrSteps <- length(Steps)
+
+ ## initialization of resulting objects:
+ ss <- vector("list", nrSteps)
+ ssRet <- vector("list", nrSteps)
+
+ withuExo <- !is.null(stateEq at uExofct)
+ withwExo <- !is.null(obsEq at wExofct)
+ withdots.prop <- (length(dots.propagated)>0)
+
+ ### noch herauszufinden: woher findet man raus,
+ ## ob der prepstep/predstep/corrstep control/Diagnostic hat..
+
+ withcontrol.smooth <- ##!is.null(stateEq at uExofct)
+ withDiagnostic.smooth <- ##!!is.null(stateEq at uExofct)
+
+ smoothret0 <- initSSPredOrFiltRet(Model at pdim, tY,
+ Model at pdim, tT, Model at qdim, tY,
+ withuExo, withwExo, withdots.prop,
+ withcontrol.smooth, withDiagnostic.smooth,
+ smooth = TRUE)
+
+ for(iStep in 1:nrStep) smoothRet[[iStep]] <- smoothret0
+
+ ### Schleife ueber Zeitpunkte t mit entsprechendem Time-Management ###
+
+ iy <- 0
+ for(ix in loopIndex){
+ ## correction:
+ if(inX[ix]){ ## have an observation available
+ iy <- iy + 1
+ for (iStep in 1:nrStep) {
+ ss[[iStep]] <- Steps[[iStep]]@smoothStep(i=iy, t=tt[ix], Obs=Obs,
+ PredOrFilt=ps[[iStep]],
+ smoothEq=smoothEq, ...)
+ ssRet[[iStep]] <- updateSSPredOrFilt(csRet[[iStep]],
+ cs[[iStep]], iy, smooth = TRUE)
+ }
+ }else{
+ cs <- ps
+ }
+ }
+ zusammenhängen von list(initRet, prepRet, psRet, csRet) ## SSOutput ist jetzt
+ eine Liste mit 4 Elementen; jedes einzelne der 4 ist wieder Liste
+ und enthält zB ideal robust -> unterschied zu Output.dia
+
+ }
\ No newline at end of file
Modified: branches/robKalman_2012/pkg/robKalman/R/updateinitSSPredOrFilterRet.R
===================================================================
--- branches/robKalman_2012/pkg/robKalman/R/updateinitSSPredOrFilterRet.R 2014-04-16 16:34:59 UTC (rev 75)
+++ branches/robKalman_2012/pkg/robKalman/R/updateinitSSPredOrFilterRet.R 2015-01-16 15:59:10 UTC (rev 76)
@@ -1,8 +1,9 @@
-updateSSPredOrFiltRet <- function(old, new, i){
+updateSSPredOrFiltRet <- function(old, new, i, smooth=FALSE){
### later: type checking
old at values[,i] <- new at values
if(!is.null(new at call)) old at call[[i]] <- new at call
old at variance[,,i] <- new at variance
+ if(smooth) old at lag1variance[,,i] <- new at lag1variance
if(!is.null(new at uExo)) old at uExo[,i] <- new at uExo
if(!is.null(new at wExo)) old at wExo[,i] <- new at wExo
if(!is.null(new at dots.propagated)) old at dots.propagated[[i]] <- new at dots.propagated
@@ -12,7 +13,7 @@
}
initSSPredOrFiltRet <- function(rdim, trdim, pdim, tdim, qdim, tydim, withuExo, withwExo, withdots.prop,
- withcontrol, withDiagnosticFilter){
+ withcontrol, withDiagnosticFilter, smooth=FALSE){
v <- matrix(NA,rdim,trdim)
vm <- array(NA,dim=c(rdim,rdim,trdim))
uExo <- if(withuExo) matrix(NA,pdim,tdim) else NULL
@@ -20,7 +21,9 @@
dots.prop <- if(withdots.prop) vector("list", trdim) else NULL
control <- if(withcontrol) vector("list", trdim) else NULL
DiagnosticFilter <- if(DiagnosticFilter) vector("list",trdim) else NULL
- new("SSPredOrFiltRet", value = v, variance = vm, uExo = uExo, wExo = wExo,
+ l1vm <- if(smooth) array(NA,dim=c(rdim,rdim,trdim-1) else NULL
+ new("SSPredOrFiltRet", value = v, variance = vm, lag1variance = l1vm,
+ uExo = uExo, wExo = wExo,
dot.propagated = dot.prop, control = control,
DiagnosticFilter = DiagnosticFilter)
}
\ No newline at end of file
More information about the Robkalman-commits
mailing list