Linear Kalman Filter Example

Introduction

This notebook is designed to demonstrate how to use the StateSpace.jl package to execute the Kalman filter for a linear State Space model. The example that has been used here closely follows the one given on "Greg Czerniak's Website". Namely the voltage example on this page.

For those of you that do not need/want the explanation of the model and the code, you can skip right to the end of this notebook where the entire section of code required to run this example is given.

The Problem

The problem that we will consider here is that of predicting the true value of the voltage given some noisy measurements from a voltmeter. We also are armed with the knowledge that the true voltage remains constant.

Process Model

We will define our state transistion as

$$x_i = x_{i-1} + v_i$$

where $x_i$ is the current state, i.e. the current voltage, $x_{i-1}$ is the previous state (voltage) and $v_i$ is the process noise error.

Observation Model

We will assume that our voltmeter gives us the voltage directly. This means the observation equation is

$$y_i = x_i + w_i$$

where $y_i$ is the current observation, i.e. measured voltage, and $w_i$ is the observation error.

Setting up the problem

First we'll import the modules that we need



In [6]:

    
using StateSpace
using Distributions
using DataFrames
using Gadfly
using Colors

Generate noisy observations

In this section we need to generate the noisy observations. Let's first define a true voltage value.



In [7]:

    
true_voltage = 1.25









    Out[7]:





1.25

Now we need to set the noise level for the observations. We'll do this by setting the variance.



In [8]:

    
measurement_noise_variance = 0.1









    Out[8]:





0.1

Now we can simulate a set of noisy measurements. We'll do this by adding zero mean Gaussian noise with the variance that we just specified above.



In [9]:

    
number_of_observations = 60 
observations = randn(number_of_observations) * sqrt(measurement_noise_variance) + true_voltage
observations = observations'









    Out[9]:





1x60 Array{Float64,2}:
 0.867216  0.94188  0.519849  0.699067  …  0.93253  1.55718  1.38088  1.10902

Note that we had to transpose the observations matrix because by default a 1-dimensional array in Julia is a column array. In StateSpace.jl a single column is considered as a single observation (for the case where we have multiple elements that are considered a single observation). This means that the second dimension is interpreted as separate observations.

Let's plot the observations to see what we are dealing with.



In [10]:

    
plt = plot(x=1:60, y=observations, Geom.point)









    Out[10]:

Define Kalman Filter Parameters

We can define the process and observation parameters according to the model defined above. Since the process doesn't change, the process parameter is equal to 1. Also since the voltage is observed directly, the observation parameter is also equal to 1. Since we are very sure that the voltage is constant we can set the process variance to a small value. We've also already set the measurement variance previously. Therefore we have the following parameters:



In [11]:

    
process_parameter = 1.0
process_variance = 0.00001
observation_parameter = 1.0
observation_variance = measurement_noise_variance









    Out[11]:





0.1

We can now create the instance of the Linear State Space Model



In [12]:

    
linSSM = LinearGaussianSSM(process_parameter, process_variance, observation_parameter, observation_variance)









    Out[12]:





StateSpace.LinearGaussianSSM{Float64}(1x1 Array{Float64,2}:
 1.0,1x1 Array{Float64,2}:
 1.0e-5,1x1 Array{Float64,2}:
 1.0,1x1 Array{Float64,2}:
 0.1)

Initial Guess

Let's create an initial guess for the voltage. We'll assume that we're not too sure and we'll make a guess that differs from the true value with some variance



In [13]:

    
initial_guess = MvNormal([3.0], [1.0])









    Out[13]:





DiagNormal(
dim: 1
μ: [3.0]
Σ: 1x1 Array{Float64,2}:
 1.0
)

Perform Kalman Filter Algorithm

Now we have all of the parameters:

noisy observations
process (transition) and observation (emission) model paramaters
initial guess of state

We can use the Kalman Filter to predict the true underlying state (voltage).



In [14]:

    
filtered_state = filter(linSSM, observations, initial_guess)









    



SmoothedState{Float64}





    Out[14]:












    



61 estimates of 1-D process from 1-D observations
Log-likelihood: 89.37290609952812

And there you have it. You have just used the StateSpace.jl package to obtain a filtered estimate of the state of the State Space system.

Plotting the result

There are several plotting packages available to this and you can find out about them here and pick your favourite one. We will demonstrate how to plot the results using Gadfly.

First we extract the filtered state along with their corresponding $2\sigma$ values (the give the area that we would expect the true value to lie with 95% confidence).



In [15]:

    
x_data = 1:number_of_observations
state_array = Vector{Float64}(number_of_observations+1)
confidence_array = Vector{Float64}(number_of_observations+1)
for i in 1:number_of_observations+1
    current_state = filtered_state.state[i]
    state_array[i] = current_state.μ[1]
    if i != 1
        confidence_array[i] = 2*sqrt(current_state.Σ.mat[1])
    else
        confidence_array[i] = 2*sqrt(current_state.Σ.diag[1])
    end
end

Next we will create a dataframe. This is simply so the syntax is simple for plotting the ribbon digram which will represent the state along with the confidence interval



In [16]:

    
df_fs = DataFrame(
    x = [0;x_data],
    y = state_array,
    ymin = state_array - confidence_array,
    ymax = state_array + confidence_array,
    f = "Filtered values"
    )









    Out[16]:




x y ymin ymax f
1 0 3.0 1.0 5.0 Filtered values
2 1 1.0611037286205813 0.45808076536638553 1.664126691874777 Filtered values
3 2 1.0043272434764778 0.5678787864387553 1.4407757005142003 Filtered values
4 3 0.8480157101346791 0.4887724921158011 1.207258928153557 Filtered values
5 4 0.8116731388775696 0.4992675786216292 1.12407869913351 Filtered values
6 5 0.807189150166465 0.5270451691959738 1.087333131136956 Filtered values
7 6 0.8999378498054265 0.6437422910175581 1.1561334085932948 Filtered values
8 7 0.989622311169391 0.7521068521282737 1.2271377702105084 Filtered values
9 8 1.0483116723177255 0.8258898085494348 1.2707335360860161 Filtered values
10 9 1.0446087784073468 0.8347087112635011 1.2545088455511924 Filtered values
11 10 1.0786302684277618 0.8793335246866631 1.2779270121688604 Filtered values
12 11 1.0884449727380594 0.898275278613649 1.2786146668624698 Filtered values
13 12 1.0982623786198673 0.9160548695418606 1.280469887697874 Filtered values
14 13 1.0841097536282016 0.9089260617901457 1.2592934454662577 Filtered values
15 14 1.0828040878336194 0.9138750565756308 1.251733119091608 Filtered values
16 15 1.0649439537758276 0.9016297312794417 1.2282581762722136 Filtered values
17 16 1.0605109946230231 0.902272468502058 1.2187495207439882 Filtered values
18 17 1.0773102444580036 0.9236881113747883 1.2309323775412189 Filtered values
19 18 1.088700606810123 0.9392997223316368 1.2381014912886092 Filtered values
20 19 1.0931892689236813 0.9476667356764505 1.238711802170912 Filtered values
21 20 1.0901126635818 0.9481686211482534 1.2320567060153464 Filtered values
22 21 1.1071350757266631 0.9685054791374939 1.2457646723158324 Filtered values
23 22 1.1189577989626351 0.9834086849427635 1.2545069129825068 Filtered values
24 23 1.1208734006302368 0.9881962833201385 1.253550517940335 Filtered values
25 24 1.1075417754310246 0.97754991100588 1.2375336398561694 Filtered values
26 25 1.1256060295414745 0.9981313573451419 1.253080701737807 Filtered values
27 26 1.1470517101934734 1.0219423264016363 1.2721610939853105 Filtered values
28 27 1.1545402778628855 1.0316583325610842 1.2774222231646868 Filtered values
29 28 1.152720888493037 1.0319408225312914 1.2735009544547828 Filtered values
30 29 1.1516722757129403 1.0328793325488863 1.2704652188769943 Filtered values
&vellip &vellip &vellip &vellip &vellip &vellip

Next we will create separate colours for the observations and the true value. Here we will generate the default distinguishable colors used by Gadfly.



In [17]:

    
n = 3 #We will require 3 different colors
getColors = distinguishable_colors(n, Color[LCHab(70, 60, 240)],
                                   transform=c -> deuteranopic(c, 0.5),
                                   lchoices=Float64[65, 70, 75, 80],
                                   cchoices=Float64[0, 50, 60, 70],
                                   hchoices=linspace(0, 330, 24))









    Out[17]:

Finally we will plot the results



In [18]:

    
filtered_state_plot = plot(
    layer(x=x_data, y=filtered_state.observations, Geom.point, Theme(default_color=getColors[2])),
    layer(x=[0;x_data], y=ones(number_of_observations+1)*true_voltage, Geom.line, Theme(default_color=getColors[3])),
    layer(df_fs, x=:x, y=:y, ymin=:ymin, ymax=:ymax, Geom.line, Geom.ribbon),
    Guide.xlabel("Measurement Number"), Guide.ylabel("Voltage (Volts)"),
    Guide.manual_color_key("Colour Key",["Filtered Estimate", "Measurements","True Value "],[getColors[1],getColors[2],getColors[3]]),
    Guide.title("Linear Kalman Filter Example")
    )
display(filtered_state_plot)

Smoothing

We can also perform smoothing on the filtered estimates to get a better estimate of the true underlying state (as they say, hindsight is 20/20). This can be performed once you have the filtered estimate like so:



In [19]:

    
smoothed_state = smooth(linSSM, filtered_state)









    



SmoothedState{Float64}





    Out[19]:












    



60 estimates of 1-D process from 1-D observations
Log-likelihood: 120.15074208037329

Now we can the results in the same way as above to visualize the smoothed result.



In [20]:

    
state_array = Vector{Float64}(number_of_observations)
confidence_array = Vector{Float64}(number_of_observations)
for i in x_data
    current_state = smoothed_state.state[i]
    state_array[i] = current_state.μ[1]
    confidence_array[i] = 2*sqrt(current_state.Σ.mat[1])
end
df_ss = DataFrame(
    x = x_data,
    y = state_array,
    ymin = state_array - confidence_array,
    ymax = state_array + confidence_array,
    f = "Filtered values"
    )

smoothed_state_plot = plot(
    layer(x=x_data, y=smoothed_state.observations, Geom.point, Theme(default_color=getColors[2], line_width=3px)),
    layer(x=x_data, y=ones(number_of_observations)*true_voltage, Geom.line, Theme(default_color=getColors[3], line_width=3px)),
    layer(df_ss, x=:x, y=:y, ymin=:ymin, ymax=:ymax, Geom.line, Geom.ribbon, Theme(line_width=3px)),
    Guide.xlabel("Measurement Number"), Guide.ylabel("Voltage (Volts)"),
    Guide.manual_color_key("Colour Key",["Smoothed Estimate", "Measurements","True Value "],[getColors[1],getColors[2],getColors[3]]),
    Guide.title("Linear Kalman Smoother Example"), Theme(major_label_font_size=20px, key_label_font_size=20px)
    )
display(smoothed_state_plot)

For those that would like to copy and paste the code. Here is the code in full:



In [21]:

    
using StateSpace
using Distributions
using DataFrames
using Gadfly
using Colors

#Generate noisy observations
true_voltage = 1.25
measurement_noise_variance = 0.1
number_of_observations = 60 
observations = randn(number_of_observations) * sqrt(measurement_noise_variance) + true_voltage
observations = observations'

#Define Filter parameters
process_parameter = 1.0
process_variance = 0.00001
observation_parameter = 1.0
observation_variance = measurement_noise_variance
linSSM = LinearGaussianSSM(process_parameter, process_variance, observation_parameter, observation_variance)

initial_guess = MvNormal([3.0], [1.0]) #Give initial guess parameters

filtered_state = filter(linSSM, observations, initial_guess) #Perform Kalman Filter
smoothed_state = smooth(linSSM, filtered_state) #Perform Kalman Smoother

#Plot the filter results
x_data = 1:number_of_observations
state_array = Vector{Float64}(number_of_observations+1)
confidence_array = Vector{Float64}(number_of_observations+1)
for i in 1:number_of_observations+1
    current_state = filtered_state.state[i]
    state_array[i] = current_state.μ[1]
    if i != 1
        confidence_array[i] = 2*sqrt(current_state.Σ.mat[1])
    else
        confidence_array[i] = 2*sqrt(current_state.Σ.diag[1])
    end
end
df_fs = DataFrame(
    x = [0;x_data],
    y = state_array,
    ymin = state_array - confidence_array,
    ymax = state_array + confidence_array,
    f = "Filtered values"
    )

n = 3
getColors = distinguishable_colors(n, Color[LCHab(70, 60, 240)],
                                   transform=c -> deuteranopic(c, 0.5),
                                   lchoices=Float64[65, 70, 75, 80],
                                   cchoices=Float64[0, 50, 60, 70],
                                   hchoices=linspace(0, 330, 24))
filtered_state_plot = plot(
    layer(x=x_data, y=filtered_state.observations, Geom.point, Theme(default_color=getColors[2])),
    layer(x=[0;x_data], y=ones(number_of_observations+1)*true_voltage, Geom.line, Theme(default_color=getColors[3])),
    layer(df_fs, x=:x, y=:y, ymin=:ymin, ymax=:ymax, Geom.line, Geom.ribbon),
    Guide.xlabel("Measurement Number"), Guide.ylabel("Voltage (Volts)"),
    Guide.manual_color_key("Colour Key",["Filtered Estimate", "Measurements","True Value "],[getColors[1],getColors[2],getColors[3]]),
    Guide.title("Linear Kalman Filter Example")
    )
display(filtered_state_plot)

#Plot the smoothed results
state_array = Vector{Float64}(number_of_observations)
confidence_array = Vector{Float64}(number_of_observations)
for i in x_data
    current_state = smoothed_state.state[i]
    state_array[i] = current_state.μ[1]
    confidence_array[i] = 2*sqrt(current_state.Σ.mat[1])
end
df_ss = DataFrame(
    x = x_data,
    y = state_array,
    ymin = state_array - confidence_array,
    ymax = state_array + confidence_array,
    f = "Filtered values"
    )

smoothed_state_plot = plot(
    layer(x=x_data, y=smoothed_state.observations, Geom.point, Theme(default_color=getColors[2], line_width=3px)),
    layer(x=x_data, y=ones(number_of_observations)*true_voltage, Geom.line, Theme(default_color=getColors[3], line_width=3px)),
    layer(df_ss, x=:x, y=:y, ymin=:ymin, ymax=:ymax, Geom.line, Geom.ribbon, Theme(line_width=3px)),
    Guide.xlabel("Measurement Number"), Guide.ylabel("Voltage (Volts)"),
    Guide.manual_color_key("Colour Key",["Smoothed Estimate", "Measurements","True Value "],[getColors[1],getColors[2],getColors[3]]),
    Guide.title("Linear Kalman Smoother Example"), Theme(major_label_font_size=20px, key_label_font_size=20px)
    )
display(smoothed_state_plot)

	x	y	ymin	ymax	f
1	0	3.0	1.0	5.0	Filtered values
2	1	1.0611037286205813	0.45808076536638553	1.664126691874777	Filtered values
3	2	1.0043272434764778	0.5678787864387553	1.4407757005142003	Filtered values
4	3	0.8480157101346791	0.4887724921158011	1.207258928153557	Filtered values
5	4	0.8116731388775696	0.4992675786216292	1.12407869913351	Filtered values
6	5	0.807189150166465	0.5270451691959738	1.087333131136956	Filtered values
7	6	0.8999378498054265	0.6437422910175581	1.1561334085932948	Filtered values
8	7	0.989622311169391	0.7521068521282737	1.2271377702105084	Filtered values
9	8	1.0483116723177255	0.8258898085494348	1.2707335360860161	Filtered values
10	9	1.0446087784073468	0.8347087112635011	1.2545088455511924	Filtered values
11	10	1.0786302684277618	0.8793335246866631	1.2779270121688604	Filtered values
12	11	1.0884449727380594	0.898275278613649	1.2786146668624698	Filtered values
13	12	1.0982623786198673	0.9160548695418606	1.280469887697874	Filtered values
14	13	1.0841097536282016	0.9089260617901457	1.2592934454662577	Filtered values
15	14	1.0828040878336194	0.9138750565756308	1.251733119091608	Filtered values
16	15	1.0649439537758276	0.9016297312794417	1.2282581762722136	Filtered values
17	16	1.0605109946230231	0.902272468502058	1.2187495207439882	Filtered values
18	17	1.0773102444580036	0.9236881113747883	1.2309323775412189	Filtered values
19	18	1.088700606810123	0.9392997223316368	1.2381014912886092	Filtered values
20	19	1.0931892689236813	0.9476667356764505	1.238711802170912	Filtered values
21	20	1.0901126635818	0.9481686211482534	1.2320567060153464	Filtered values
22	21	1.1071350757266631	0.9685054791374939	1.2457646723158324	Filtered values
23	22	1.1189577989626351	0.9834086849427635	1.2545069129825068	Filtered values
24	23	1.1208734006302368	0.9881962833201385	1.253550517940335	Filtered values
25	24	1.1075417754310246	0.97754991100588	1.2375336398561694	Filtered values
26	25	1.1256060295414745	0.9981313573451419	1.253080701737807	Filtered values
27	26	1.1470517101934734	1.0219423264016363	1.2721610939853105	Filtered values
28	27	1.1545402778628855	1.0316583325610842	1.2774222231646868	Filtered values
29	28	1.152720888493037	1.0319408225312914	1.2735009544547828	Filtered values
30	29	1.1516722757129403	1.0328793325488863	1.2704652188769943	Filtered values
&vellip	&vellip	&vellip	&vellip	&vellip	&vellip