Observed variables

[dfabianus]

Hey @baggepinnen ✋

Is it possible to obtain observed variables from the estimation that are not part of the state vector?

My usecase is the following (state transition function and output function):

function FLUMO_discrete(x,u,p,t; tS=0.001)
    V, I, N, A, D, k_offset1, k_offset2  = x
    n, P0, F0, beta1, b_n, a_n, b_a, a_a = p
    c_D = D ./ V
    k_n = maximum([a_n .* (1 .+ abs(c_D)) .^ b_n .+ k_offset1, 0])
    k_a = maximum([a_a .* (1 .+ abs(c_D)) .^ b_a .+ k_offset2, 0])
    dI = -k_n .* I .- k_a .* I .^ n
    dN = k_n .* I
    dA = k_a .* I .^ n
    Vn = x[1] .+ u[3] #Volume with discrete pulse events
    In = x[2] .+ u[2] .+ dI .* tS #Intermediates with discrete pulse events
    Nn = x[3] .+ dN .* tS #Native protein
    An = x[4] .+ dA .* tS #Aggregated protein
    Dn = x[5] .+ u[1] #Denaturant with discrete pulse events
    return [Vn, In, Nn, An, Dn, k_offset1, k_offset2]
end

function FLUMO_discrete_measurement(x,u,p,t)
    V, I, N, _, D, k_offset1,k_offset2  = x
    n, P0, F0, beta1, b_n, a_n, b_a, a_a = p
    c_D = D ./ V
    k_n = maximum([a_n .* (1 .+ abs(c_D)) .^ b_n .+ k_offset1, 0])
    k_a = maximum([a_a .* (1 .+ abs(c_D)) .^ b_a .+ k_offset2, 0])
    dI = -k_n .* (I./V) .- k_a .* (I./V) .^ n
    dAEW = dI ./ beta1
    F = ((I+N)./V)*F0/(P0./V)
    return [F, dAEW]
end

I would be interested in visualizing the trajectories of k_n and k_a which are algebraically computed inside the model functions but are not part of the state vector. k_offset1 and k_offset2 however are part of the state vector which are estimating the error of k_n and k_a. I assume that "observed variables" are not possible out of the box in the numeric formulation.

Therefore I think I have the following options:

• Recompute k_n and k_a after state estimation from the obtained estimates but this seems a little inefficient
• Introduce k_n and k_a into the state vector with zero noise. I assume that is the best approach.

Do you have an advice for the best practice in such a case?

Thanks a lot!

[baggepinnen]

I would go the route recomputing the k_n, k_a variables here. The expressions for those are very cheap to compute compared to the computations that take place inside the filter. An UKF or EKF will perform matrix factorizations that scale as $\mathcal{O}(n_x^3)$, so adding trivial state variables to the dynamics is very expensive.

A side note,

maximum([a_n .* (1 .+ abs(c_D)) .^ b_n .+ k_offset1, 0])

is more efficiently computed using

max(a_n * (1 + abs(c_D)) ^ b_n + k_offset1, 0)

This avoids allocating an array.

I see that you use broadcasting in a lot of places where I'd expect scalars to appear, is there a reason for this? For example, why c_D = D ./ V instead of c_D = D / V? Using broadcasting when it's not needed increases compile times and can also fail to catch mistakes.

[dfabianus]

Yes, when I think about it, it makes sense that more states increase the complexity even more.

Thank you also for the other recommendations ;)