Missing Observations Example

Introduction

This example is a really basic/unrealistic example that was made up to demonstrate the use of missing observations with the StateSpace.jl package.

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.

NOTE: An explanation of the affects of the different missing observation implementations is given just prior to the final chunk of code.

The Problem

In this example we have three (spherical) balloons that deflate over time. However we only record (noisy) measurements of the radius of only one of the balloons over time. We consider 3 cases for the unobserved balloon:

Our best guess is the prior prediction
We have some information on the measurement error of the balloons (even the unobserved ones)
We know that there is some correlation between two of the balloons.

The Process Model

We assume the process is that the balloon's radius decreases by 10% of it's previous radius with each observation. Written mathematically, this is

$$ \mathbf{x}_i = 0.9 \times \mathbf{x}_{i-1} + \mathbf{v}_i, $$

Where $\mathbf{x}_i$ represents the state at time step $i$ and $\mathbf{v}_i$ is the process noise error at time step $i$. The state here is the vector containing the radius of each balloon.

The Observation Model

We assume that we observe the radius of the the balloon directly, which means the observation equation is

$$\mathbf{y}_i = \mathbf{x}_i + \mathbf{w}_i$$

where $\mathbf{y}_i$ is the current observation, i.e. the radius, and $\mathbf{w}_i$ is the observation error.

Setting up the problem

First we'll import the required modules for performing the filtering algorithm and visualising the results



In [1]:

    
using StateSpace
using Distributions
using Gadfly
using DataFrames
using Colors

Set Guess of Initial State

Here we begin with a guess of the initial state of the system i.e. what the initial radii of the balloons are. Because we don't know any better, we'll set the variances to 1.



In [6]:

    
initialState = [30.0, 15.0, 15.0]
initialCov = eye(length(initialState))
initialGuess = MvNormal(initialState, initialCov)









    Out[6]:





FullNormal(
dim: 3
μ: [30.0,15.0,15.0]
Σ: 3x3 Array{Float64,2}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0
)

Generate Noisy Observations

Here we'll generate the noisy observations of each balloon's radius. Notice that there are only measurements for the first balloon - We don't observe the other two balloons. To represent missing observations we use NaN values in their place. So here is the code to generate the noisy observations along with the missing ones.



In [7]:

    
numObs = 30
numBalloons = length(initialState)
observationVariance = 2.0
observations = Matrix{Float64}(numBalloons,numObs)
trueValue = Vector{Float64}(numObs)
trueValue[1] = 0.9 * initialState[1]
observations[1,1] = trueValue[1] + randn()*sqrt(observationVariance)
observations[2:end,1] = [NaN, NaN]
for i in 2:numObs
    trueValue[i] = 0.9*trueValue[i-1]
    observations[1,i] = trueValue[i] + randn()*sqrt(observationVariance)
    observations[2:end,i] = [NaN, NaN]
end

Describe Kalman Filter parameters

We now need to define the process and observation model according to the model described above.



In [8]:

    
#Define the process matrix and covariance matrix.
processMatrix = 0.9 * eye(initialCov)
processCov = eye(initialCov)

observationMatrix = eye(processMatrix) #Define the observation matrix









    Out[8]:





3x3 Array{Float64,2}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

For this example we're going to explicitly describe the relationship of the measurements bewteen the balloons with a correlation matrix (note that a correlation matrix is a normalized covariance matrix with values between -1 and 1 corresponding to perfect anticorrelation and perfect correlation respectively). However the Kalman filter works with covariance matrices. This means that we need to define a function to convert bewteen a correlation matrix and a covariance matrix. So we'll define this function now.



In [9]:

    
function corr2covar{T}(correlationMatrix::Matrix{T}, sigmaVector::Vector{T})
    N = length(sigmaVector)
    covarianceMatrix = zeros(correlationMatrix)
    for j in 1:N
        for i in 1:N
            covarianceMatrix[i, j] = correlationMatrix[i,j] * sigmaVector[i] * sigmaVector[j]
        end
    end
    return covarianceMatrix
end









    Out[9]:





corr2covar (generic function with 1 method)

Now we can define a correlation matrix and the standard deviations of the balloon measurement errors which can be converted into a covariance matrix. In this matrix, we assume the balloon measurements are uncorrelated. (For case 3 we'll change this to ensure two of the balloons have correlated measurement errors).



In [10]:

    
#Correlation Matrix
observationCorr = [1.0 0.0 0.0;
                   0.0 1.0 0.0;
                   0.0 0.0 1.0]
observationΣ = [1.0, 3.0, 3.0] #measurement error standard deviations
observationCov = corr2covar(observationCorr, observationΣ) #Conversion to covariance matrix









    Out[10]:





3x3 Array{Float64,2}:
 1.0  0.0  0.0
 0.0  9.0  0.0
 0.0  0.0  9.0

Now we have defined all of the parameters of our model we can create the linear state space model object



In [11]:

    
linSSM = LinearGaussianSSM(processMatrix, processCov, observationMatrix, observationCov)









    Out[11]:





StateSpace.LinearGaussianSSM{Float64}(3x3 Array{Float64,2}:
 0.9  0.0  0.0
 0.0  0.9  0.0
 0.0  0.0  0.9,3x3 Array{Float64,2}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0,3x3 Array{Float64,2}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0,3x3 Array{Float64,2}:
 1.0  0.0  0.0
 0.0  9.0  0.0
 0.0  0.0  9.0)

Case 1: Our best guess is the prior prediction

If we're completely uncertain of the measurement variances then we update the state with our predicted state. The error in doing this will propagte and we'll have an estimate with a relatively large uncertainty. The syntax for this case is the same as the usual syntax for the package since this is the default behaviour (and the safest, although not necessarily optimal).



In [12]:

    
filtState_c1 = filter(linSSM, observations, initialGuess)









    



SmoothedState{Float64}





    Out[12]:












    



31 estimates of 3-D process from 3-D observations
Log-likelihood: -150.78794786777019

Case 2: We have some information on the measurement error of the balloons (even the unobserved ones)

We have some information on the measurement error. If we think that we have some information about the measurement errors then we can take advantage of that by using it on our predicted observersation (i.e. using our measurement errors on our predicted observations that are generated by applying the observation model to the state that is predicted from the process model). The syntax for this requires an extra parameter for the filter method - 'true'. This tells the filter that we are happy to use our measurement error estimates.



In [13]:

    
filtState_c2 = filter(linSSM, observations, initialGuess, true)









    



SmoothedState{Float64}





    Out[13]:












    



31 estimates of 3-D process from 3-D observations
Log-likelihood: -289.3515218741301

Case 3: We know that there is some correlation between two of the balloons.

For this case we assume that there is some correlation between the measurement errors. So we need to redefine the observation covariance matrix and recreate the linear statespace model object.



In [15]:

    
observationCorr = [1.0 0.6 0.0;
                   0.6 1.0 0.0;
                   0.0 0.0 1.0]
observationCov = corr2covar(observationCorr, observationΣ)
linSSM2 = LinearGaussianSSM(processMatrix, processCov, observationMatrix, observationCov)









    Out[15]:





StateSpace.LinearGaussianSSM{Float64}(3x3 Array{Float64,2}:
 0.9  0.0  0.0
 0.0  0.9  0.0
 0.0  0.0  0.9,3x3 Array{Float64,2}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0,3x3 Array{Float64,2}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0,3x3 Array{Float64,2}:
 1.0  1.8  0.0
 1.8  9.0  0.0
 0.0  0.0  9.0)

Now we can perform the Kalman filter with the correlated observation errors. Notice that the syntax is exactly the same as before.



In [16]:

    
filtState_c3 = filter(linSSM2, observations, initialGuess, true)









    



SmoothedState{Float64}





    Out[16]:












    



31 estimates of 3-D process from 3-D observations
Log-likelihood: -282.04471199474114

Plot the results

Here we'll plot the results of the fitering runs. The first chunk of code is concerned with extracting the appropriate values and creating the dataframe objects for each balloon in each of the cases



In [17]:

    
#Balloon radius and 95% confidence for each balloon - case 1
b1_c1 = Vector{Float64}(numObs+1)
b2_c1 = Vector{Float64}(numObs+1)
b3_c1 = Vector{Float64}(numObs+1)
b1Sig_c1 = Vector{Float64}(numObs+1)
b2Sig_c1 = Vector{Float64}(numObs+1)
b3Sig_c1 = Vector{Float64}(numObs+1)

#Balloon radius and 95% confidence for each balloon - case 2
b1_c2 = Vector{Float64}(numObs+1)
b2_c2 = Vector{Float64}(numObs+1)
b3_c2 = Vector{Float64}(numObs+1)
b1Sig_c2 = Vector{Float64}(numObs+1)
b2Sig_c2 = Vector{Float64}(numObs+1)
b3Sig_c2 = Vector{Float64}(numObs+1)

#Balloon radius and 95% confidence for each balloon - case 3
b1_c3 = Vector{Float64}(numObs+1)
b2_c3 = Vector{Float64}(numObs+1)
b3_c3 = Vector{Float64}(numObs+1)
b1Sig_c3 = Vector{Float64}(numObs+1)
b2Sig_c3 = Vector{Float64}(numObs+1)
b3Sig_c3 = Vector{Float64}(numObs+1)

#Extract state estimates from each balloon for each case
for i in 1:numObs+1
    currentState1 = filtState_c1.state[i]
    b1_c1[i] = currentState1.μ[1]
    b2_c1[i] = currentState1.μ[2]
    b3_c1[i] = currentState1.μ[3]
    b1Sig_c1[i] = 2*sqrt(currentState1.Σ.mat[1,1])
    b2Sig_c1[i] = 2*sqrt(currentState1.Σ.mat[2,2])
    b3Sig_c1[i] = 2*sqrt(currentState1.Σ.mat[3,3])
    
    currentState2 = filtState_c2.state[i]
    b1_c2[i] = currentState2.μ[1]
    b2_c2[i] = currentState2.μ[2]
    b3_c2[i] = currentState2.μ[3]
    b1Sig_c2[i] = 2*sqrt(currentState2.Σ.mat[1,1])
    b2Sig_c2[i] = 2*sqrt(currentState2.Σ.mat[2,2])
    b3Sig_c2[i] = 2*sqrt(currentState2.Σ.mat[3,3])
    
    currentState3 = filtState_c3.state[i]
    b1_c3[i] = currentState3.μ[1]
    b2_c3[i] = currentState3.μ[2]
    b3_c3[i] = currentState3.μ[3]
    b1Sig_c3[i] = 2*sqrt(currentState3.Σ.mat[1,1])
    b2Sig_c3[i] = 2*sqrt(currentState3.Σ.mat[2,2])
    b3Sig_c3[i] = 2*sqrt(currentState3.Σ.mat[3,3])
end

#Create data frames for each balloon for each case
df_b1_c1 = DataFrame(
    x = 1:numObs+1,
    y = b1_c1,
    ymin = b1_c1 - b1Sig_c1,
    ymax = b1_c1 + b1Sig_c1,
    f = "Filtered values"
    )

df_b2_c1 = DataFrame(
    x = 1:numObs+1,
    y = b2_c1,
    ymin = b2_c1 - b2Sig_c1,
    ymax = b2_c1 + b2Sig_c1,
    f = "Filtered values"
    )

df_b3_c1 = DataFrame(
    x = 1:numObs+1,
    y = b3_c1,
    ymin = b3_c1 - b3Sig_c1,
    ymax = b3_c1 + b3Sig_c1,
    f = "Filtered values"
    )

df_b1_c2 = DataFrame(
    x = 1:numObs+1,
    y = b1_c2,
    ymin = b1_c2 - b1Sig_c2,
    ymax = b1_c2 + b1Sig_c2,
    f = "Filtered values"
    )

df_b2_c2 = DataFrame(
    x = 1:numObs+1,
    y = b2_c2,
    ymin = b2_c2 - b2Sig_c2,
    ymax = b2_c2 + b2Sig_c2,
    f = "Filtered values"
    )

df_b3_c2 = DataFrame(
    x = 1:numObs+1,
    y = b3_c2,
    ymin = b3_c2 - b3Sig_c2,
    ymax = b3_c2 + b3Sig_c2,
    f = "Filtered values"
    )

df_b1_c3 = DataFrame(
    x = 1:numObs+1,
    y = b1_c3,
    ymin = b1_c3 - b1Sig_c3,
    ymax = b1_c3 + b1Sig_c3,
    f = "Filtered values"
    )

df_b2_c3 = DataFrame(
    x = 1:numObs+1,
    y = b2_c3,
    ymin = b2_c3 - b2Sig_c3,
    ymax = b2_c3 + b2Sig_c3,
    f = "Filtered values"
    )

df_b3_c3 = DataFrame(
    x = 1:numObs+1,
    y = b3_c3,
    ymin = b3_c3 - b3Sig_c3,
    ymax = b3_c3 + b3Sig_c3,
    f = "Filtered values"
    )









    Out[17]:




x y ymin ymax f
1 1 15.0 13.0 17.0 Filtered values
2 2 13.5 11.04485192271299 15.95514807728701 Filtered values
3 3 12.15 9.480807286921358 14.819192713078642 Filtered values
4 4 10.935 8.16280281751531 13.707197182484691 Filtered values
5 5 9.8415 7.019591014724654 12.663408985275346 Filtered values
6 6 8.85735 6.01145258370574 11.70324741629426 Filtered values
7 7 7.971615000000001 5.114146538037939 10.829083461962062 Filtered values
8 8 7.174453500000001 4.3114051075634645 10.037501892436538 Filtered values
9 9 6.457008150000001 3.5912693229701245 9.322746977029878 Filtered values
10 10 5.811307335000001 2.9442713757028165 8.678343294297186 Filtered values
11 11 5.230176601500001 2.3625152818344324 8.09783792116557 Filtered values
12 12 4.707158941350001 1.8391961345584806 7.575121748141521 Filtered values
13 13 4.236443047215001 1.368334894295177 7.104551200134825 Filtered values
14 14 3.8127987424935013 0.9446205188884051 6.680976966098598 Filtered values
15 15 3.4315188682441513 0.5633068639585841 6.299730872529718 Filtered values
16 16 3.0883669814197363 0.2201386916751633 5.956595271164309 Filtered values
17 17 2.779530283277763 -0.0887058575831281 5.647766424138654 Filtered values
18 18 2.5015772549499866 -0.3666626708833509 5.369817180783324 Filtered values
19 19 2.251419529454988 -0.6168222210889116 5.119661279998888 Filtered values
20 20 2.0262775765094894 -0.8419650537154415 4.89452020673442 Filtered values
21 21 1.8236498188585404 -1.0445932354548522 4.691892873171933 Filtered values
22 22 1.6412848369726865 -1.226958421790973 4.509528095736346 Filtered values
23 23 1.477156353275418 -1.3910870040523744 4.3453997106032105 Filtered values
24 24 1.3294407179478762 -1.5388026868970375 4.19768412279279 Filtered values
25 25 1.1964966461530886 -1.6717467815995177 4.0647400739056945 Filtered values
26 26 1.0768469815377797 -1.7913964572584749 3.9450904203340342 Filtered values
27 27 0.9691622833840018 -1.8990811607363232 3.837405727504327 Filtered values
28 28 0.8722460550456016 -1.9959973916414229 3.740489501732626 Filtered values
29 29 0.7850214495410415 -2.0832219983833715 3.653264897465455 Filtered values
30 30 0.7065193045869373 -2.1617241439340136 3.574762753107888 Filtered values
&vellip &vellip &vellip &vellip &vellip &vellip

Now we obtain the 3 colours required for each balloon



In [18]:

    
#Get the correct colours for the plots
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))









    Out[18]:



In [22]:

    
plt_c1 = plot(
    layer(x=2:numObs+1, y=observations[1,:], Geom.point, Theme(default_color=getColors[1])),
    layer(df_b1_c1, x=:x, y=:y, ymin=:ymin, ymax=:ymax, Geom.line, Geom.ribbon, Theme(default_color=getColors[1])),
    layer(df_b2_c1, x=:x, y=:y, ymin=:ymin, ymax=:ymax, Geom.line, Geom.ribbon, Theme(default_color=getColors[2])),
    layer(df_b3_c1, x=:x, y=:y, ymin=:ymin, ymax=:ymax, Geom.line, Geom.ribbon, Theme(default_color=getColors[3])),
    Guide.xlabel("Measurement Number"), Guide.ylabel("Balloon Radius (cm)"),
    Guide.manual_color_key("Colour Key",["Balloon 1", "Balloon 2","Balloon 3"],[getColors[1],getColors[2],getColors[3]]),
    Guide.title("Missing Observations Example - Case 1")
    )
display(plt_c1)



In [20]:

    
plt_c2 = plot(
    layer(x=2:numObs+1, y=observations[1,:], Geom.point, Theme(default_color=getColors[1])),
    layer(df_b1_c2, x=:x, y=:y, ymin=:ymin, ymax=:ymax, Geom.line, Geom.ribbon, Theme(default_color=getColors[1])),
    layer(df_b2_c2, x=:x, y=:y, ymin=:ymin, ymax=:ymax, Geom.line, Geom.ribbon, Theme(default_color=getColors[2])),
    layer(df_b3_c2, x=:x, y=:y, ymin=:ymin, ymax=:ymax, Geom.line, Geom.ribbon, Theme(default_color=getColors[3])),
    Guide.xlabel("Measurement Number"), Guide.ylabel("Balloon Radius (cm)"),
    Guide.manual_color_key("Colour Key",["Balloon 1", "Balloon 2","Balloon 3"],[getColors[1],getColors[2],getColors[3]]),
    Guide.title("Missing Observations Example - Case 2")
    )
display(plt_c2)



In [24]:

    
plt_c3 = plot(
    layer(x=2:numObs+1, y=observations[1,:], Geom.point, Theme(default_color=getColors[1])),
    layer(df_b1_c3, x=:x, y=:y, ymin=:ymin, ymax=:ymax, Geom.line, Geom.ribbon, Theme(default_color=getColors[1])),
    layer(df_b2_c3, x=:x, y=:y, ymin=:ymin, ymax=:ymax, Geom.line, Geom.ribbon, Theme(default_color=getColors[2])),
    layer(df_b3_c3, x=:x, y=:y, ymin=:ymin, ymax=:ymax, Geom.line, Geom.ribbon, Theme(default_color=getColors[3])),
    Guide.xlabel("Measurement Number"), Guide.ylabel("Balloon Radius (cm)"),
    Guide.manual_color_key("Colour Key",["Balloon 1", "Balloon 2","Balloon 3"],[getColors[1],getColors[2],getColors[3]]),
    Guide.title("Missing Observations Example - Case 3")
    )
display(plt_c3)

Explanation of the Results

Case 1

In the first plot (case 1) we have the case where rather than make any inference about the state of the system using any measurements we may have, we simply just assume that the new state is exactly what we predicted - that is the case for all three balloons. This is made more obvious by looking at the blue curve which doesn't use any information from the measurements. Notice in this case that the errors on the estimate are propagated and so the shaded region (which represents the 95% confidence interval) gets larger as time goes on. This is the default behaviour of this package. This is because it's the safest option if you have no information at all about missing measurements. As highlighted in this example, all elements in the state vector are assumed to follow the predicted value, even if there are partial measurements at a given time interval i.e. in this case, even though we observed measurements for balloon 1, we didn't use those measurements because the other balloons had missing values.

Case 2

In this case (the second plot) we notice that the blue line actually seems to use information from the measurements. This is because we have told the filter to use any measurements that we are given. If there are missing values then these observations are estimated. (so make sure you know what you are doing if you use this option). Since the balloon measurement errors are uncorrelated, the radius for balloons 2 and 3 are exactly the same (again). Their confidence region is smaller than in case 1 (i.e. we're more sure about the radius), simply because of the size of the measurement error standard deviation we gave earlier. If these values were increased for balloon 2 and/or 3, then we would see an increase in the interval. So I reiterate, if you use this option, you have to make sure you know what you're doing. You need an idea about the measurement error that you expect to have given your prediction.

Case 3

This case uses exactly the same algorithm as case two - i.e. any of the actual observations are used as normal and the missing observations are estimated. But in this case the measurement errors of balloon 1 and balloon 2 were correlated. This means that we use this information to update the state of balloon 2. In the third plot you can see that balloon 2 deviates from what would usually be predicted (i.e. it deviates from the estimate of balloon 3). You can play about with the correlation matrix defined earlier to investigate how the correlation affects the state estimate.

For those that just want the code, without the explanation, you have come to the right place. Here it is:



In [1]:

    
using StateSpace
using Distributions
using Gadfly
using DataFrames
using Colors

#Set guess of the initial state
initialState = [30.0, 15.0, 15.0]
initialCov = eye(length(initialState))
initialGuess = MvNormal(initialState, initialCov)


#Generate noisy Observations
numObs = 30
numBalloons = length(initialState)
observationVariance = 2.0
observations = Matrix{Float64}(numBalloons,numObs)
trueValue = Vector{Float64}(numObs)
trueValue[1] = 0.9 * initialState[1]
observations[1,1] = trueValue[1] + randn()*sqrt(observationVariance)
observations[2:end,1] = [NaN, NaN]
for i in 2:numObs
    trueValue[i] = 0.9*trueValue[i-1]
    observations[1,i] = trueValue[i] + randn()*sqrt(observationVariance)
    observations[2:end,i] = [NaN, NaN]
end

#Describe Kalman Filter parameters
processMatrix = 0.9 * eye(initialCov)
processCov = eye(initialCov)

#Define function to convert a correlation matrix to a covariance matrix 
function corr2covar{T}(correlationMatrix::Matrix{T}, sigmaVector::Vector{T})
    N = length(sigmaVector)
    covarianceMatrix = zeros(correlationMatrix)
    for j in 1:N
        for i in 1:N
            covarianceMatrix[i, j] = correlationMatrix[i,j] * sigmaVector[i] * sigmaVector[j]
        end
    end
    return covarianceMatrix
end
observationMatrix = eye(processMatrix) #Define observation matrix

#define the correlation matrix
observationCorr = [1.0 0.0 0.0;
                   0.0 1.0 0.0;
                   0.0 0.0 1.0]
#Now we define a vector containing standard deviation values
observationΣ = [1.0, 3.0, 3.0]
observationCov = corr2covar(observationCorr, observationΣ) #convert correlation matrix to covariance matrix

#Create the linear StateSpace model.
linSSM = LinearGaussianSSM(processMatrix, processCov,
                           observationMatrix, observationCov)


#Section: Perform Kalman Filter for cases 1 and 2
#Case 1: Our best guess is the prior prediction
filtState_c1 = filter(linSSM, observations, initialGuess)

#Case 2: We have some information on the measurement error
filtState_c2 = filter(linSSM, observations, initialGuess, true)

#Section: Perform Kalman Filter for case 3
#Case 3: We know that there is some correlation between balloons 1 and 2.
observationCorr = [1.0 0.6 0.0;
                   0.6 1.0 0.0;
                   0.0 0.0 1.0]
observationCov = corr2covar(observationCorr, observationΣ) #convert correlation matrix to covariance matrix
#Create the linear StateSpace model with new observation covariance matrix.
linSSM2 = LinearGaussianSSM(processMatrix, processCov,
                           observationMatrix, observationCov)

filtState_c3 = filter(linSSM2, observations, initialGuess, true)

#Plot Filtered results
#Balloon radius and 95% confidence for each balloon - case 1
b1_c1 = Vector{Float64}(numObs+1)
b2_c1 = Vector{Float64}(numObs+1)
b3_c1 = Vector{Float64}(numObs+1)
b1Sig_c1 = Vector{Float64}(numObs+1)
b2Sig_c1 = Vector{Float64}(numObs+1)
b3Sig_c1 = Vector{Float64}(numObs+1)

#Balloon radius and 95% confidence for each balloon - case 2
b1_c2 = Vector{Float64}(numObs+1)
b2_c2 = Vector{Float64}(numObs+1)
b3_c2 = Vector{Float64}(numObs+1)
b1Sig_c2 = Vector{Float64}(numObs+1)
b2Sig_c2 = Vector{Float64}(numObs+1)
b3Sig_c2 = Vector{Float64}(numObs+1)

#Balloon radius and 95% confidence for each balloon - case 3
b1_c3 = Vector{Float64}(numObs+1)
b2_c3 = Vector{Float64}(numObs+1)
b3_c3 = Vector{Float64}(numObs+1)
b1Sig_c3 = Vector{Float64}(numObs+1)
b2Sig_c3 = Vector{Float64}(numObs+1)
b3Sig_c3 = Vector{Float64}(numObs+1)

#Extract state estimates from each balloon for each case
for i in 1:numObs+1
    currentState1 = filtState_c1.state[i]
    b1_c1[i] = currentState1.μ[1]
    b2_c1[i] = currentState1.μ[2]
    b3_c1[i] = currentState1.μ[3]
    b1Sig_c1[i] = 2*sqrt(currentState1.Σ.mat[1,1])
    b2Sig_c1[i] = 2*sqrt(currentState1.Σ.mat[2,2])
    b3Sig_c1[i] = 2*sqrt(currentState1.Σ.mat[3,3])
    
    currentState2 = filtState_c2.state[i]
    b1_c2[i] = currentState2.μ[1]
    b2_c2[i] = currentState2.μ[2]
    b3_c2[i] = currentState2.μ[3]
    b1Sig_c2[i] = 2*sqrt(currentState2.Σ.mat[1,1])
    b2Sig_c2[i] = 2*sqrt(currentState2.Σ.mat[2,2])
    b3Sig_c2[i] = 2*sqrt(currentState2.Σ.mat[3,3])
    
    currentState3 = filtState_c3.state[i]
    b1_c3[i] = currentState3.μ[1]
    b2_c3[i] = currentState3.μ[2]
    b3_c3[i] = currentState3.μ[3]
    b1Sig_c3[i] = 2*sqrt(currentState3.Σ.mat[1,1])
    b2Sig_c3[i] = 2*sqrt(currentState3.Σ.mat[2,2])
    b3Sig_c3[i] = 2*sqrt(currentState3.Σ.mat[3,3])
end

#Create data frames for each balloon for each case
df_b1_c1 = DataFrame(
    x = 1:numObs+1,
    y = b1_c1,
    ymin = b1_c1 - b1Sig_c1,
    ymax = b1_c1 + b1Sig_c1,
    f = "Filtered values"
    )

df_b2_c1 = DataFrame(
    x = 1:numObs+1,
    y = b2_c1,
    ymin = b2_c1 - b2Sig_c1,
    ymax = b2_c1 + b2Sig_c1,
    f = "Filtered values"
    )

df_b3_c1 = DataFrame(
    x = 1:numObs+1,
    y = b3_c1,
    ymin = b3_c1 - b3Sig_c1,
    ymax = b3_c1 + b3Sig_c1,
    f = "Filtered values"
    )

df_b1_c2 = DataFrame(
    x = 1:numObs+1,
    y = b1_c2,
    ymin = b1_c2 - b1Sig_c2,
    ymax = b1_c2 + b1Sig_c2,
    f = "Filtered values"
    )

df_b2_c2 = DataFrame(
    x = 1:numObs+1,
    y = b2_c2,
    ymin = b2_c2 - b2Sig_c2,
    ymax = b2_c2 + b2Sig_c2,
    f = "Filtered values"
    )

df_b3_c2 = DataFrame(
    x = 1:numObs+1,
    y = b3_c2,
    ymin = b3_c2 - b3Sig_c2,
    ymax = b3_c2 + b3Sig_c2,
    f = "Filtered values"
    )

df_b1_c3 = DataFrame(
    x = 1:numObs+1,
    y = b1_c3,
    ymin = b1_c3 - b1Sig_c3,
    ymax = b1_c3 + b1Sig_c3,
    f = "Filtered values"
    )

df_b2_c3 = DataFrame(
    x = 1:numObs+1,
    y = b2_c3,
    ymin = b2_c3 - b2Sig_c3,
    ymax = b2_c3 + b2Sig_c3,
    f = "Filtered values"
    )

df_b3_c3 = DataFrame(
    x = 1:numObs+1,
    y = b3_c3,
    ymin = b3_c3 - b3Sig_c3,
    ymax = b3_c3 + b3Sig_c3,
    f = "Filtered values"
    )

#Get the correct colours for the plots
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))

#Plot the balloon radius estimates for each case
plt_c1 = plot(
    layer(x=2:numObs+1, y=observations[1,:], Geom.point, Theme(default_color=getColors[1])),
    layer(df_b1_c1, x=:x, y=:y, ymin=:ymin, ymax=:ymax, Geom.line, Geom.ribbon, Theme(default_color=getColors[1])),
    layer(df_b2_c1, x=:x, y=:y, ymin=:ymin, ymax=:ymax, Geom.line, Geom.ribbon, Theme(default_color=getColors[2])),
    layer(df_b3_c1, x=:x, y=:y, ymin=:ymin, ymax=:ymax, Geom.line, Geom.ribbon, Theme(default_color=getColors[3])),
    Guide.xlabel("Measurement Number"), Guide.ylabel("Balloon Radius (cm)"),
    Guide.manual_color_key("Colour Key",["Balloon 1", "Balloon 2","Balloon 3"],[getColors[1],getColors[2],getColors[3]]),
    Guide.title("Missing Observations Example - Case 1")
    )
display(plt_c1)

plt_c2 = plot(
    layer(x=2:numObs+1, y=observations[1,:], Geom.point, Theme(default_color=getColors[1])),
    layer(df_b1_c2, x=:x, y=:y, ymin=:ymin, ymax=:ymax, Geom.line, Geom.ribbon, Theme(default_color=getColors[1])),
    layer(df_b2_c2, x=:x, y=:y, ymin=:ymin, ymax=:ymax, Geom.line, Geom.ribbon, Theme(default_color=getColors[2])),
    layer(df_b3_c2, x=:x, y=:y, ymin=:ymin, ymax=:ymax, Geom.line, Geom.ribbon, Theme(default_color=getColors[3])),
    Guide.xlabel("Measurement Number"), Guide.ylabel("Balloon Radius (cm)"),
    Guide.manual_color_key("Colour Key",["Balloon 1", "Balloon 2","Balloon 3"],[getColors[1],getColors[2],getColors[3]]),
    Guide.title("Missing Observations Example - Case 2")
    )
display(plt_c2)

plt_c3 = plot(
    layer(x=2:numObs+1, y=observations[1,:], Geom.point, Theme(default_color=getColors[1])),
    layer(df_b1_c3, x=:x, y=:y, ymin=:ymin, ymax=:ymax, Geom.line, Geom.ribbon, Theme(default_color=getColors[1])),
    layer(df_b2_c3, x=:x, y=:y, ymin=:ymin, ymax=:ymax, Geom.line, Geom.ribbon, Theme(default_color=getColors[2])),
    layer(df_b3_c3, x=:x, y=:y, ymin=:ymin, ymax=:ymax, Geom.line, Geom.ribbon, Theme(default_color=getColors[3])),
    Guide.xlabel("Measurement Number"), Guide.ylabel("Balloon Radius (cm)"),
    Guide.manual_color_key("Colour Key",["Balloon 1", "Balloon 2","Balloon 3"],[getColors[1],getColors[2],getColors[3]]),
    Guide.title("Missing Observations Example - Case 3")
    )
display(plt_c3)

	x	y	ymin	ymax	f
1	1	15.0	13.0	17.0	Filtered values
2	2	13.5	11.04485192271299	15.95514807728701	Filtered values
3	3	12.15	9.480807286921358	14.819192713078642	Filtered values
4	4	10.935	8.16280281751531	13.707197182484691	Filtered values
5	5	9.8415	7.019591014724654	12.663408985275346	Filtered values
6	6	8.85735	6.01145258370574	11.70324741629426	Filtered values
7	7	7.971615000000001	5.114146538037939	10.829083461962062	Filtered values
8	8	7.174453500000001	4.3114051075634645	10.037501892436538	Filtered values
9	9	6.457008150000001	3.5912693229701245	9.322746977029878	Filtered values
10	10	5.811307335000001	2.9442713757028165	8.678343294297186	Filtered values
11	11	5.230176601500001	2.3625152818344324	8.09783792116557	Filtered values
12	12	4.707158941350001	1.8391961345584806	7.575121748141521	Filtered values
13	13	4.236443047215001	1.368334894295177	7.104551200134825	Filtered values
14	14	3.8127987424935013	0.9446205188884051	6.680976966098598	Filtered values
15	15	3.4315188682441513	0.5633068639585841	6.299730872529718	Filtered values
16	16	3.0883669814197363	0.2201386916751633	5.956595271164309	Filtered values
17	17	2.779530283277763	-0.0887058575831281	5.647766424138654	Filtered values
18	18	2.5015772549499866	-0.3666626708833509	5.369817180783324	Filtered values
19	19	2.251419529454988	-0.6168222210889116	5.119661279998888	Filtered values
20	20	2.0262775765094894	-0.8419650537154415	4.89452020673442	Filtered values
21	21	1.8236498188585404	-1.0445932354548522	4.691892873171933	Filtered values
22	22	1.6412848369726865	-1.226958421790973	4.509528095736346	Filtered values
23	23	1.477156353275418	-1.3910870040523744	4.3453997106032105	Filtered values
24	24	1.3294407179478762	-1.5388026868970375	4.19768412279279	Filtered values
25	25	1.1964966461530886	-1.6717467815995177	4.0647400739056945	Filtered values
26	26	1.0768469815377797	-1.7913964572584749	3.9450904203340342	Filtered values
27	27	0.9691622833840018	-1.8990811607363232	3.837405727504327	Filtered values
28	28	0.8722460550456016	-1.9959973916414229	3.740489501732626	Filtered values
29	29	0.7850214495410415	-2.0832219983833715	3.653264897465455	Filtered values
30	30	0.7065193045869373	-2.1617241439340136	3.574762753107888	Filtered values
&vellip	&vellip	&vellip	&vellip	&vellip	&vellip