Real-time data reconciliation and gross error detection

[ValentinKaisermayer]

Is it possible to do real-time data reconciliation and gross error detection with RxInfer?

E.g. I have y1, y2 and y3 as noisy measurements and know that y1 + y2 = y3. Or If there is some redundancy in the measurements but one measures with non-zero mean error to detect it.

[wouterwln]

What you are trying to do should be possible in RxInfer! From what I understand, your model looks something like this:

@model function error_correction(y)
    x[1] ~ NormalMeanVariance(y[1], 1)
    x[2] ~ NormalMeanVariance(y[2], 1)
    x[3] ~ NormalMeanVariance(y[3], 1)
    x[3] := x[1] + x[2]
end

Here we have some noisy observations y, and there are "true" latents x which are related through the equation you have listed. In this model, inference is quite straightforward, and with infer(model = error_correction(), data = (y = [1, 2, 3],)) you will get posterior estimates over your latents x. However, you might want to infer the "noise" introduced by your measurement device, so let's go and put a prior on that:

@model function error_correction(y, n_1_shape, n_1_rate, n_2_shape, n_2_rate, n_3_shape, n_3_rate)
    n_1 ~ GammaShapeRate(n_1_shape, n_1_rate)
    x[1] ~ NormalMeanPrecision(y[1], n_1)
    n_2 ~ GammaShapeRate(n_2_shape, n_2_rate)
    x[2] ~ NormalMeanPrecision(y[2], n_2)
    n_3 ~ GammaShapeRate(n_3_shape, n_3_rate)
    x[3] ~ NormalMeanPrecision(y[3], n_3)
    x[3] := x[1] + x[2]
end

I've added a lot of arguments to this model that govern the prior distribution over the noise distribution. That is because I want to move towards the real-time part of your question. I want to incorporate the example here , such that this model can perform inference on an infinite stream of data. For this we also need to define how the prior over noise is updated over time. Since we are also inferring noise of a Gaussian distribution now, we have to initialize constraints, and we have to give an initial value for the noise parameters in order to kickstart inference:

updates = @autoupdates begin
    n_1_shape = shape(q(n_1))
    n_1_rate = rate(q(n_1))
    n_2_shape = shape(q(n_2))
    n_2_rate = rate(q(n_2))
    n_3_shape = shape(q(n_3))
    n_3_rate = rate(q(n_3))
end

initialization = @initialization begin
    q(n_1) = GammaShapeRate(1, 1)
    q(n_2) = GammaShapeRate(1, 1)
    q(n_3) = GammaShapeRate(1, 1)
end

constraints = @constraints begin
    q(n_1, n_2, n_3,  x) = MeanField()
end
RxInfer.is_data(x::Vector{RxInfer.GraphVariableRef}) = RxInfer.is_data(first(x))

I snuck this last line in there to make this example work. This is a known issue in RxInfer and my fix is quite an ugly one. This should be resolved in the next version of RxInfer so you can ignore it for now. Now we are in a shape where we can generate some data:

observations = [[rand(Normal(0, 1)) for _ in 1:2] for _ in 1:100]
for observation in observations
    push!(observation, rand(Normal(sum(observation), 0.3)))
end

# create a stream from this static dataset, in your pipeline this would come from some real-time sensors, this is just to illustrate that RxInfer can handle these datastreams
static_datastream   = from(observations) |> map(NamedTuple{(:y,), Tuple{Vector{Float64}}}, (d) -> (y = d, ));

Now a call to infer will give you posterior estimates over your noise parameters and over the true value of your measurements:

infer(model = error_correction(), 
     datastream = static_datastream, 
     autoupdates = updates, 
     initialization = initialization, 
     constraints = constraints
)

Additional options and visualisation options can be found here. I hope this answers your question! If there's anything else, feel free to ask.

[ValentinKaisermayer]

Thank you, that is really helpful!
Do I need to use a hierarchical model to be able to use infinite datastreams and @autoupdate? Let's say I know the variance in my measurements. Would @autoupdate be empty in this case? I known that I could use a weighted least squares approach in that case.

Can you give an example of an implementation when I have a dynamical model? E.g.
x[k+1] = x[k] + T * (y1[k] - y2[k]), where this is just a simple integrator, with T being the sample time of the data. I would have measurements of x, y1 and y2.

[wouterwln]

I think your model, from a probabilistic sense, is a bit ill-defined; You define x[k+1] recursively using x[k], but in the next timestep you say that you observe x[k+1] as well, which would mean that you assign every x twice. Maybe x is latent and you want to combine beliefs from both a "predictive" distribution and a measurement? This would turn the model into a regular state filtering problem (with a slightly different likelihood but because of the message passing nature of RxInfer this wouldn't make much of a difference) which can be done with a similar approach as the Kalman Filter example in the documentation.

The main idea is that in @autoupdates you specify which information of the posterior distribution you want to carry over as priors when the next datapoint arrives (which would be your state if I interpret your model correctly). You do not have to use a hierarchical model to use infinite datastreams. In fact, @autoupdates can only supply state updates to the outermost model in your hierarchical model, so in this context it might be beneficial to fit everything in a single @model block.

[ValentinKaisermayer]

I meant that I have noisy measurements of x[k], u1[k] and u2[k]. My latent variables should fulfil the model given the noisy observations.

@model function error_correction(y,x)
    x_lat[1] ~ NormalMeanVariance(x[1],1) 
    u1_lat[1] ~ NormalMeanVariance(u1[1],1) 
    u2_lat[1] ~ NormalMeanVariance(u2[1],1) 

    for t in 2:length(y)
        u1_lat[t] ~ NormalMeanVariance(u1[t],1) 
        u2_lat[t] ~ NormalMeanVariance(u2[t],1)
        x_lat[t] := x_lat[t-1] + tS * (u1_lat[t-1] + u2_lat[t-1])
        x[t] ~ NormalMeanVariance(x_lat[t],1) 
    end
end

In terms of the KF-example, I would like to estimate the inputs u to the system.
An application would be if I have a sensor fault for u1 for some time and have missing values, I would still be able to estimate u1_lat, right?

[wouterwln]

Yes. What you can do is that you keep track of the latent states of u and x_lat. This can be fed with @autoupdates to a next iteration of the model, which will then learn from x to disentangle the two contributions. If you pass missings, indeed RxInfer will do predictions for these variables instead of passing messages backwards to the latent variables, so you could handle the missing data by estimating the forward direction of your generative model. Also, in your model you set y as an input to the model, but you only use it in the call to length. In a streaming data example, you wouldn't need this for loop at all so I think you can drop y from the list of interfaces.

[ValentinKaisermayer]

Sry, it should be u

@model function error_correction(u,x)
    x_lat[1] ~ NormalMeanVariance(x[1],1) 
    u1_lat[1] ~ NormalMeanVariance(u[1][1],1) 
    u2_lat[1] ~ NormalMeanVariance(u[2][1],1) 

    for t in 2:length(u)
        u1_lat[t] ~ NormalMeanVariance(u[1][t],1) 
        u2_lat[t] ~ NormalMeanVariance(u[2][t],1)
        x_lat[t] := x_lat[t-1] + tS * (u1_lat[t-1] + u2_lat[t-1])
        x[t] ~ NormalMeanVariance(x_lat[t],1) 
    end
end

For the autoupdate:

@model function error_correction(x_prev_mean, x_prev_var, u1_prev_mean, u1_prev_var, u2_prev_mean, u2_prev_var, u_prev, x)
    x_prev ~ Normal(mean = x_prev_mean, variance = x_prev_var)
    u1_prev ~ Normal(mean = u1_prev_mean, variance = u1_prev_var)
    u2_prev ~ Normal(mean = u2_prev_mean, variance = u2_prev_var)
    x_curent := x_prev + tS * (u1_prev + u2_prev)
    x ~ Normal(mean = x_curent , precision = 1)  # observe `x`
    # where do I observe `u`?
end

autoupdates = @autoupdates begin 
    x_prev_mean, x_prev_var = mean_var(q(x_curent))
    u1_prev_mean, u1_prev_var = mean_var(q(u1_prev))
    u2_prev_mean, u2_prev_var = mean_var(q(u2_prev))
end

[ValentinKaisermayer]

keep track of the latent states of u

Why should I keep track of u? It is not a static parameter that I want to update my beliefs of. In your first example you had the unknown variance of the measurements but these are parameters. Do I understand this correctly?

[wouterwln]

I think you observe u in this line: u1_lat[t] ~ NormalMeanVariance(u[1][t],1). Because of the reflexivity of the Gaussian distribution in its mean parameter you could also write: u[1] ~ NormalMeanVariance(u1_prev, 1), right? As for your second comment, you're right. You should probably only keep track of u1 and u2 as you have done and in your model explain how they transition between datapoints.