All observation schemes inherit from
ObservationSchemes.Observation
which has methods
ObservationSchemes.eltype
ObservationSchemes.size
ObservationSchemes.length
automatically implemented for it. The idea is to decorate each recorded data-point with such structs, and in doing so, encode the way in which it was collected.
We implemented two concrete structs that may be used for defining a single observation:
LinearGsnObs
: to encode observations of the linear transformations of the underlying process disturbed by Gaussian noiseGeneralObs
: to encode observations of non-linear transformations of the underlying process disturbed by noise.
Both admit a possibility of parameterization.
The most important observation scheme is LinearGsnObs
. It is suitable for representing observations that can be written in the following format:
where
ObservationSchemes.LinearGsnObs
Below, we list some special cases of the scheme above.
A degenerate case of the setting above is an exact observation of
so that
t, v = 1.0, [1.0, 2.0, 3.0]
obs = LinearGsnObs(t, v; full_obs=true)
!!! warning
For numerical reasons the covariance matrix of the noise Σ
should not be a zero matrix, and instead, even in the exact observation setting should be inflated by some small artificial noise
. By default Σ
is set to
Σ
explicitly, for instance to increase the level of artificial noise
:
julia using LinearAlgebra obs = LinearGsnObs(t, v; Σ=(1e-5)I, full_obs=true)
Specifying full_obs=true
is important, as it lets the compiler differentiate
between an actual, full observation with some artificial noise and a (possibly)
partial observation with very low level (or also artificial level) of noise.
As a result, other packages from JuliaDiffusionBayes know when the Markov property can be applied.
We can view a summary of the observation by calling summary
:
julia> summary(obs)
⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤
|Observation `v = Lx+ξ`, where `L` is a (3, 3)-matrix, `x` is a state of the stochastic process and `ξ`∼N(μ,Σ).
|...
|| v: [1.0, 2.0, 3.0] (observation),
|| → typeof(v): Array{Float64,1},
|| made at time 1.0.
|...
|L: [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0],
| → typeof(L): Diagonal{Float64,Array{Float64,1}}
|μ: [0.0, 0.0, 0.0],
| → typeof(μ): Array{Float64,1}
|Σ: [1.0e-11 0.0 0.0; 0.0 1.0e-11 0.0; 0.0 0.0 1.0e-11],
| → typeof(Σ): Diagonal{Float64,Array{Float64,1}}
|...
|This is an exact observation.
|...
|It does not depend on any additional parameters.
|...
|No first passage times recorded.
⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆
Notice that various defaults and type-inferences have kicked in. It was recognized that the observation does not depend on any parameters, that first-passage time setting does not apply and that the observation was not of a static type and hence regular Arrays
are used to define L
, μ
and Σ
.
This is a standard understanding of the expression in \eqref{eq:obs_scheme}. An example could be:
using StaticArrays
v = @SVector [1.0, 2.0, 3.0]
t = 2.0
L = @SMatrix [1.0 0.0 2.0 0.0; 3.0 4.0 0.0 0.0; 0.0 1.0 0.0 1.0]
Σ = SDiagonal(1.0, 1.0, 1e-11)
obs = LinearGsnObs(t, v; L = L, Σ = Σ)
for defining a three-dimensional observation v
made at time
julia> summary(obs)
⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤
|Observation `v = Lx+ξ`, where `L` is a (3, 4)-matrix, `x` is a state of the stochastic process and `ξ`∼N(μ,Σ).
|...
|| v: [1.0, 2.0, 3.0] (observation),
|| → typeof(v): SArray{Tuple{3},Float64,1,3},
|| made at time 2.0.
|...
|L: [1.0 0.0 2.0 0.0; 3.0 4.0 0.0 0.0; 0.0 1.0 0.0 1.0],
| → typeof(L): SArray{Tuple{3,4},Float64,2,12}
|μ: [0.0, 0.0, 0.0],
| → typeof(μ): SArray{Tuple{3},Float64,1,3}
|Σ: [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0e-11],
| → typeof(Σ): Diagonal{Float64,SArray{Tuple{3},Float64,1,3}}
|...
|This is NOT an exact observation.
|...
|It does not depend on any additional parameters.
|...
|No first passage times recorded.
⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆
Note that the internal containers are now set to be SVector
s (even μ
, which wasn't passed to a constructor but its type was inferred and its value set to zero). Additionally, Julia understands that this is not a full observation and hence Markov property cannot be applied.
Support for certain first-passage time settings is provided. By default LinearGsnObs
sets the information about the first-passage times to
ObservationSchemes.NoFirstPassageTimes
However, this can be changed by passing appropriately initialized
ObservationSchemes.FirstPassageTimeInfo
For instance, to indicate in the example above that the last coordinate of v
actually reaches level
t, v = 1.0, [1.0, 2.0, 3.0]
fpt = FirstPassageTimeInfo(
(3,),
(3.0,),
(true,),
(false,),
(),
)
obs = LinearGsnObs(t, v; L = L, Σ = Σ, fpt = fpt)
The last two entries in fpt
specify additional reset times. For instance,
instead of (false,)
and ()
we could set (true,)
, (-1.0)
to indicate that
the process
julia> summary(obs)
⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤
|Observation `v = Lx+ξ`, where `L` is a (3, 4)-matrix, `x` is a state of the stochastic process and `ξ`∼N(μ,Σ).
|...
|| v: [1.0, 2.0, 3.0] (observation),
|| → typeof(v): Array{Float64,1},
|| made at time 1.0.
|...
|L: [1.0 0.0 2.0 0.0; 3.0 4.0 0.0 0.0; 0.0 1.0 0.0 1.0],
| → typeof(L): SArray{Tuple{3,4},Float64,2,12}
|μ: [0.0, 0.0, 0.0],
| → typeof(μ): Array{Float64,1}
|Σ: [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0e-11],
| → typeof(Σ): Diagonal{Float64,SArray{Tuple{3},Float64,1,3}}
|...
|This is NOT an exact observation.
|...
|It does not depend on any additional parameters.
|...
|First passage times of the observation `v`:
|----------------------------------------------------------------------------
|| coordinate | level | up-crossing | extra reset | reset level |
|----------------------------------------------------------------------------
|| 3 | 3.0 | up-crossing | ✘ | ✘ |
⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆
changes appropriately to display the new summary of the first-passage time information.
!!! note This package is agnostic with respect to the algorithms that are later used on the decorated observations. Consequently it doesn't make any checks for whether the observations make sense. For instance in the package GuidedProposals.jl that deals with simulating conditioned diffusions, a support for first-passage time observations is currently extended only to diffusions where the dynamics of the coordinate whose first-passage time is observed is devoid of any Wiener noise. The onus of checking whether this or other constraints are satisfied are on the user.
!!! tip "Why do we refer to LinearGsnObs
as most important?"
In practice, all other observation schemes are handled by other packages in JuliaDiffusionBayes by approximating them with a suitable LinearGsnObs
and then correcting the resulting approximation error with the Metropolis-Hastings algorithm. Consequently, LinearGsnObs
will be a building block of any other observation scheme.
The observations can be parameterized by passing a vector of parameters θ
. Additionally, a Tag
needs to be attached that is used to differentiate at compile time between the non-parameterized observations and parameterized observations as well as among different parameterizations themselves.
For instance, to indicate in the [second example](@ref standard_example_lingsnobs) that two entries in the L
matrix are parameterized we can write:
v = @SVector [1.0, 2.0, 3.0]
t = 2.0
L = @SMatrix [1.0 0.0 -99.9 0.0; 3.0 4.0 0.0 0.0; 0.0 -99.9 0.0 1.0]
Σ = SDiagonal(1.0, 1.0, 1e-11)
obs = LinearGsnObs(t, v; L = L, Σ = Σ, θ=[2.0, 1.0], Tag=1)
We are not done yet, in this case matrix L
that we passed above is incomplete as we intend to create an actual matrix L
by combining the matrix and the parameters we've passed. To this end, we must overwrite the behaviour of the function L(⋅)
:
const OBS = ObservationSchemes
function OBS.L(o::LinearGsnObs{1})
_L = MMatrix{3,4,Float64}(o.L)
_L[1,3] = o.θ[1]
_L[3,2] = o.θ[2]
SMatrix{3,4,Float64}(_L)
end
Notice that we dispatch on the observation's tag 1
. Furthermore, when calling summarize
matrix L
is displayed correctly.
julia> summary(obs)
⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤
|Observation `v = Lx+ξ`, where `L` is a (3, 4)-matrix, `x` is a state of the stochastic process and `ξ`∼N(μ,Σ).
|...
|| ν: [1.0, 2.0, 3.0] (observation),
|| → typeof(ν): SArray{Tuple{3},Float64,1,3},
|| made at time 2.0.
|...
|L: [1.0 0.0 2.0 0.0; 3.0 4.0 0.0 0.0; 0.0 1.0 0.0 1.0],
| → typeof(L): SArray{Tuple{3,4},Float64,2,12}
|μ: [0.0, 0.0, 0.0],
| → typeof(μ): SArray{Tuple{3},Float64,1,3}
|Σ: [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0e-11],
| → typeof(Σ): Diagonal{Float64,SArray{Tuple{3},Float64,1,3}}
|...
|This is NOT an exact observation.
|...
|It depends on additional parameters, which are set to: (2.0, 1.0).
|...
|No first passage times recorded.
⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆
However, now we can change parameters to new values and the matrix L
will be updated:
OBS.update_params!(obs, [7.0, 8.0])
julia> summary(obs)
⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤
|Observation `v = Lx+ξ`, where `L` is a (3, 4)-matrix, `x` is a state of the stochastic process and `ξ`∼N(μ,Σ).
|...
|| ν: [1.0, 2.0, 3.0] (observation),
|| → typeof(ν): SArray{Tuple{3},Float64,1,3},
|| made at time 2.0.
|...
|L: [1.0 0.0 7.0 0.0; 3.0 4.0 0.0 0.0; 0.0 8.0 0.0 1.0],
| → typeof(L): SArray{Tuple{3,4},Float64,2,12}
|μ: [0.0, 0.0, 0.0],
| → typeof(μ): SArray{Tuple{3},Float64,1,3}
|Σ: [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0e-11],
| → typeof(Σ): Diagonal{Float64,SArray{Tuple{3},Float64,1,3}}
|...
|This is NOT an exact observation.
|...
|It depends on additional parameters, which are set to: (7.0, 8.0).
|...
|No first passage times recorded.
⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆
!!! note
We used number -99.9
just to remind ourselves that entries with this number need to be overwritten by parameters. There are however no rules or enforcements as to how the user deals with such entries.
Notice that trying to access L
via obs.L
will give you an incorrect result. To avoid such mistakes, always query L
, μ
,Σ
and obs
via accessors:
ObservationSchemes.L
ObservationSchemes.μ
ObservationSchemes.Σ
ObservationSchemes.ν
ObservationSchemes.obs
In principle, any observation types are supported, but this comes at a cost of having to provide some information explicitly. The main struct for specifing general observation schemes is
ObservationSchemes.GeneralObs
For this struct the required parameters are the time at which the observation was recorded and the observation itself, as well as the approximation via LinearGsnObs
. For instance, to specify the following observational scheme:
where
we may write:
using Distributions
# recording
t, v = 1.5, [1.0, 2.0]
# for the observation scheme
μ, Σ = [-1.0, 2.0], [1.0 0.0; 0.0 1.0]
dist = MvTDist(4, μ, Σ)
g(x) = view(x, 1:2).^2
# for some poor, ad-hoc Gaussian approximation
L = [2.0 0.0 0.0; 0.0 2.0 0.0]
# define observation
obs = GeneralObs(t, v, LinearGsnObs(t, v; L=L, Σ=Σ, μ=μ); dist=dist, g=g)
We can now view the summary:
julia> summary(obs)
⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤⏤
|Observation `v = g(x)+ξ`, where `g` is an operator defined by a function typeof(g), and `ξ` is a random variable given byDistributions.GenericMvTDist{Float64,PDMats.PDMat{Float64,Array{Float64,2}},Array{Float64,1}}.
|...
|| ν: [1.0, 2.0, 3.0] (observation),
|| → typeof(ν): SArray{Tuple{3},Float64,1,3},
|| made at time 2.0.
|...
|This is NOT an exact observation.
|...
|It does not depend on any additional parameters.
|...
|No first passage times recorded.
|...
|To inspect the linearized approximation to this observation scheme please type in:
| summary(<name-of-the-variable>.lin_obs)
|and hit ENTER.
⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆
!!! tip
The GeneralObs
can be decorated with first-passage time information and parameters in the same way as LinearGsnObs
can. However, by design, you cannot set full_obs
to true
for GeneralObs
.