Framework for Varying Time steps

[waTeim]

This is more of a question than an error. I have an application where I want velocity combined with position similar to you linear example. However, delta t is not a constant. There is some information about what to do in this situation in this SO question and this research paper and a MATLAB inplementation. But it's unclear to me how to apply these things using StateSpace. Can you suggest a path?

[ElOceanografo]

The simplest way to deal with a varying timestep might be to create a new data array by inserting NaNs in between your unequally spaced observations. For instance, if you've got observations xa, xb, and xc at times 1, 3, and 6, the new array would be [xa, NaN, xb, NaN, NaN, xc;].

A second method would be to use a time-varying matrix where the velocity is multiplied by the appropriate timestep. For instance, with a state vector that is [x, v] (i.e. 1-D motion), it would be something like this:

const dt = [1.2, 1.0, 4.3, 0.5,...and so on and so forth]
const v = 3.0
F(t) = [1 v * dt[i]; 0, 1]

Depending on the mechanics of your situation, you might want to make the process-noise matrix time-varying too. For instance, variance is proportional to elapsed time for random-walk processes.

[waTeim]

Hey good deal; I have a followup. I think that inserting NaN won't work because the sample arrival rate is not going to be an integer multiple of any known value beforehand. But it is possible to track the time when the data point is recorded, so approach 2 is possible, however, in the example code is the following

process_matrix = [[1.0, Δt, 0.0, 0.0] [0.0, 1.0, 0.0, 0.0] [0.0, 0.0, 1.0, Δt] [0.0, 0.0, 0.0, 1.0]]'
...
linCISMM = LinearGaussianCISSM(process_matrix, process_covariance, observation_matrix, observation_covariance, control_matrix, control_input)

I imagine process_matrix would be replaced by an F(t) analog. Does LinearGaussianCISSM accept functions rather than static matrices as a parameter? Or alternately, does the model need to be recomputed for each sample?

[ElOceanografo]

The model types are actually under development and will be changing soon, and the LinearGaussianCISSM type is going away. The /dev branch has the preview. But yes, in the new version you can pass in either a matrix F or a function F(t).

[waTeim]

Is this something that will be available, or that is currently available on dev?

[ElOceanografo]

Just merged the dev branch, so time-varying matrices are now available. When you construct your LinearGaussianSSM, just pass it a function F(t) that returns a matrix of the appropriate size. The argument t is meant to be the integer timestep.

[waTeim]

Can do! and updating now and reworking. I have a further complication, though, or maybe it makes it more simple? The algorithm is online that is the samples are acquired in real-time, and I'm obtaining a new prediction for each new sample. Does F(t) need to still take a parameter t? How are the functions predict and update affected?

[ElOceanografo]

Just made another change to time-varying matrices (0fee542); they now take a Real value for the parameter t, which now represents the real continuous time, not the integer timestep. If you want your process or observation matrices to vary in time, you define a function F(t) that returns the appropriate matrix and pass that to the model constructor. If they don't vary, just pass in matrices. Under the hood, they will be converted to functions which accept an argument t, but don't do anything with it--they just return the constant matrix you gave. This isn't documented well yet, but the code is at lines 45-52 here. The predict and update methods have an optional argument t (defaults to 0.0) which is passed on to the matrix-generating functions defined in the model.

If I understand your application, you've got new measurements arriving at irregular intervals, and want to update the state (a position and velocity?) as each one comes in? If that's the case, you can now slightly hack this t argument to accomplish what you're looking for. Just define your F(t) function so that it treats t like Δt:

F(Δt) = [1.0, Δt, 0.0, 0.0; 
             0.0, 1.0, 0.0, 0.0;
             0.0, 0.0, 1.0, Δt;
             0.0, 0.0, 0.0, 1.0]
# define other functions here
mod = LinearGaussianSSM(F, V, G, W)

Then, as measurements come in, you can keep track of the time between them and pass that to predict and update.

y, this_time = get_new_measurement() # or wherever the measurement and time come from
Δt = this_time - prev_time
state_estimate = update(mod, state_estimate, y, t=Δt)
prev_time = this_time

[waTeim]

I was using this before the change and it worked, and I'm using it now and it still works. I don't have a lot of time to go into it right now. Previously I used thunks and set some global variable reference in the closure (a giant hack). That did do the trick (and still does apparently). I do want to get to your expected approach if possible. Here's a link to the code for reference.

for future, allow a option argument (a tuple) to be supplied as the arguments to the callback.