Smoother for gaussian_dlm_obs_lpdf

General
1

Quick question

Has anyone written a smoother for extracting smoothed states when using gaussian_dlm_obs_lpdf?
I see in the math library there is gaussian_dlm_obs_rng but It’s not in the documentation.

Many Thanks

3 Likes
2

I think @tadej was the last one to touch the code, so maybe he knows more and has time to answer?

3

I don’t remember working on this function, but even if I did, it was just performance tuning. So I don’t have an answer.

4

hi @betascoo8 did you manage to find a smoother for gaussian_dlm_obs_lpdf? I am working on a problem where I might need to do a Kalman smoothing. Any luck with using gaussian_dlm_obs_rng ?

5

Not really, I ended up updating / testing the work of Jeff.
I managed to get it to run, but performance wasn’t great for my application and I put it to one side. I’m sure i’m just hitting the limits of my coding skills more than anything.

For my application the core Stan function it’s useful because I have missing values I need the kalman filter to estimate anyway.

They’re rough and probably full of bugs, but hopefully the attached files and below code are helpful.

require(dlm)
require(rstan)
require(data.table)
example = dlmRandom(m = 4, p = 3, nobs = 1000)

  # // matrix [p, n] y = observations matrix
  # // matrix [p, m] Z = design or measurement matrix (F in dlm)
  # // matrix [m, m] T = transition matrix (G in dlm)
  # // matrix [p, p] H = observation covariance matrix (V in dlm)
  # // matrix [r, r] Q = state covariance (W in dlm)
  # // vector [m] a_1 = initial state mean (m0 in dlm)
  # // matrix [m, m] P_1 = initial state covariance (C0 in dlm)

y_ = as.matrix(example$y) # [p, n]
Z_ = example$mod$FF
T_ = example$mod$GG
H_ = example$mod$V
Q_ = example$mod$W
a_1 = example$mod$m0
P_1 = example$mod$C0
n_ = nrow(y_) # number of observations
p_ = ncol(y_) # number of observed states
m_ = ncol(T_) # number of states
r_ = nrow(Q_) # number of state covariance (Jeff uses q)

mod = stan("ssm_lpdf_smooth.stan",
           data = list(n = n_, m = m_, p = p_, r = r_, y = y_list,
                       d = list(d_), Z = list(Z_), H = list(H_),
                       c = list(c_), T = list(T_), R = list(R_), Q = list(Q_),
                       a1 = a_1, P1 = P_1), algorithm = "Fixed_param", chains = 1, iter = 1)
draws = extract(mod, "alpha")
out = data.table(melt(draws))

ssm_lpdf.stan (5.9 KB)
ssm_lpdf_miss.stan (971 Bytes)
ssm_lpdf_smooth.stan (986 Bytes)
ssm.stan (76.6 KB)

2 Likes
6

thanks @betascoo8. I think I am going to try using other APIs for time series state space models as my application is very simple. I might revisit some of Jeff’s and your codes if I need to do something more sophisticated that other packages don’t allow for.