MARSS icon indicating copy to clipboard operation
MARSS copied to clipboard

MARSSkfas when V0 not 0

Open eeholmes opened this issue 5 years ago • 4 comments

Version 3.11.5 Info 9/27/2021

This is completed but before closing completely do the following.

  • [ ] Check the code in last comment.
  • [ ] Search for V0 not 0 in the documentation and add a section on V0!=0 if that is missing.
  • [ ] Look at MARSSinfo() and check entry for V0!=0

-Bug. Only made it happen for Z as column vector with V0 non-zero. Guessing it is in MARSSkfss() and seems to be in the smoother part.

dat <- t(harborSealWA)
dat <- dat[2:4, ] # remove the year row
# fit a model with 1 hidden state and 3 observation time series
kemfit <- MARSS(dat, model = list(
  Z=matrix(1,3,1),
  x0=matrix(0),
  V0=matrix(10000)
))
MARSSkfas(kemfit)$VtT[,,1:2]
MARSSkfss(kemfit)$VtT[,,1:2]

Gives

 MARSSkfas(kemfit)$VtT[,,1:2]
[1] 0.006583545 0.012637968
 MARSSkfss(kemfit)$VtT[,,1:2]
[1] 0.00658371 0.01263808

With Z identity, ok

kemfit <- MARSS(dat, model = list(
  x0=matrix(0,3,1),
  V0=diag(10000,3)
))
 diag(MARSSkfas(kemfit)$VtT[,,2])
[1] 0.04039229 0.01485954 0.01070933
 diag(MARSSkfss(kemfit)$VtT[,,2])
[1] 0.04039229 0.01485954 0.01070933

eeholmes avatar Aug 17 '20 21:08 eeholmes

MARSSkfas(kemfit, return.lag.one=FALSE)$VtT[,,1:2] is the same so not in the stack part of code.

Happens with tinitx=1 or tinitx=0

eeholmes avatar Aug 17 '20 21:08 eeholmes

Not a bug per se but numerical inaccuracy in the classic Kalman filter. When Z is a column vector, an inversion of Z %*% Vtt1 %*% t(Z) + R is done. If Vtt1 is much larger than R then that is ill-conditioned. With V0 = 10000, Vtt1 at t=1 is much larger than R and this is ill-conditioned.

MARSSkfss is not used in the fitting functions. Addressed by changing the default kf function in the residuals to use MARSSkfas.

eeholmes avatar Aug 17 '20 22:08 eeholmes

Good enough. Part of the numerical problems with the classic filter/smoother. Avoid making the V0 so bad by not trying to make it diffuse.

eeholmes avatar Aug 18 '20 02:08 eeholmes

Try this for Kt[,,1] which worked for the V0T problem

  • solve(t(Vtt1[,,1]), B%*%V00)

eeholmes avatar Aug 18 '20 03:08 eeholmes