baggepinnen / lowlevelparticlefilters.jl Goto Github PK

State estimation, smoothing and parameter estimation using Kalman and particle filters.

Home Page: https://baggepinnen.github.io/LowLevelParticleFilters.jl/stable

License: Other

Julia 100.00%

particle-filter kalman-filter estimation state-estimation control-systems dynamical-systems sequential-monte-carlo bayesian-inference system-identification monte-carlo-methods unscented-kalman-filter extended-kalman-filter parameter-estimation prediction-error-method state-estimation-algorithms state-estimation-filters data-assimilation virtual-sensors controls virtual-sensing

lowlevelparticlefilters.jl's Introduction

LowLevelParticleFilters

This is a library for state estimation, smoothing and parameter estimation.

Estimator Types

We provide a number of filter types

ParticleFilter: This filter is simple to use and assumes that both dynamics noise and measurement noise are additive.
AuxiliaryParticleFilter: This filter is identical to ParticleFilter, but uses a slightly different proposal mechanism for new particles.
AdvancedParticleFilter: This filter gives you more flexibility, at the expense of having to define a few more functions.
KalmanFilter. A standard Kalman filter. Has the same features as the particle filters, but is restricted to linear dynamics (possibly time varying) and Gaussian noise.
SqKalmanFilter. A standard Kalman filter on square-root form (slightly slower but more numerically stable with ill-conditioned covariance).
ExtendedKalmanFilter: For nonlinear systems, the EKF runs a regular Kalman filter on linearized dynamics. Uses ForwardDiff.jl for linearization. The noise model must be Gaussian.
UnscentedKalmanFilter: The Unscented kalman filter often performs slightly better than the Extended Kalman filter but may be slightly more computationally expensive. The UKF handles nonlinear dynamics and measurement model, but still requires an additive Gaussian noise model.
DAEUnscentedKalmanFilter: An Unscented Kalman filter for differential-algebraic systems (DAE).

Documentation

lowlevelparticlefilters.jl's People

Contributors

Stargazers

Watchers

Forkers

yakir12 akbek balinus astupidbear laranea hcarlsso wagberg rafeemusabbir mfkiwl jankap ven-k n0wis schminin danscr justcodebase

lowlevelparticlefilters.jl's Issues

More flexible noise in UKF

It would be nice to have a version of the UKF that can handle dynamics on the form $f(x, u, w, p, t)$ where both $w$ and uncertainty in $p$ are treated by creating sigmapoints for $[x, w, p]$ instead of just adding $R_1$ after the propagation of sigmapoints based on $x$ only.

This would allow modeling non-additive noise as well as uncertain parameters. It's unclear if using this to model uncertain parameters would offer any improvement over adding the uncertain parameter as a state to be estimated, but it would be useful for the non-additive noise.

Implementation of backwards filtering, forward guiding

Reference implementation of backwards filtering, forward guiding with https://arxiv.org/abs/1712.03807
https://gist.github.com/mschauer/6f8e0f57ce0b15eef9f86390d0631df4
@mschauer

Ensemble Kalman filters?

Hello!

Thanks for your work. It's a really nice package. 😄

I was wondering if there was a way to do/create an Ensemble Kalman Filter with the current state of the library?

Observed variables

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!

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

UKF

Use sigmapoints from MCM
store seedpoints, proceed as kalman filter

Particle filter with second order Markov model

Is it possible to model the dynamics of a particle filter as
$$x_t = a_1 x_{t-1} + a_2 x_{t-2} + noise,$$
where $a_i$ are weighting factors fulfilling $\sum a_i = 1$?
I assume it should be possible to have the state be a vector $[x_{t-1}, x_{t-2}]$, but this creates quite some redundancy.
Is there a way to access the previous state from within the dynamics function?

notation + documentation

I'm having trouble figuring out what the inputs correspond to in the Kalman filter example. Can you TeX out the system or provide a reference to the text whose notation you followed? The code provides no additional help as it looks like so:

"""
KalmanFilter(A,B,C,D,R1,R2,d0=MvNormal(R1))
"""
function KalmanFilter(A,B,C,D,R1,R2,d0=MvNormal(R1))
    cR1 = cond(R1)
    cR2 = cond(R2)
    (cond(cR1) > 1e8 || cond(cR2) > 1e8) && @warn("Covariance matrices are poorly conditioned")
    
    KalmanFilter(A,B,C,D,R1,R2,MvNormal(R2), d0, Vector(d0.μ), Matrix(d0.Σ), Ref(1))
end

Dependency Deprecated - Yeppp

I am trying to install the library, but I get errors regarding Yeppp. I believe Yeppp has been deprecated and archived.
The error I receive suggest that I should try building Yeppp, but when I try to download it, the server returns error 404 Not Found.

May be you can get rid of the Yeppp dependency?

I am using Julia 1.5.0 with Intel MKL build.

I get the following error:

julia> using LowLevelParticleFilters
[ Info: Precompiling LowLevelParticleFilters [d9d29d28-c116-5dba-9239-57a5fe23875b]
ERROR: LoadError: Yeppp not properly installed. Please run Pkg.build("Yeppp")
Stacktrace:
[1] error(::String) at .\error.jl:33
[2] top-level scope at C:\Users\Sannan.julia\packages\Yeppp\LvmJH\src\Yeppp.jl:11
[3] include(::Function, ::Module, ::String) at .\Base.jl:380
[4] include(::Module, ::String) at .\Base.jl:368
[5] top-level scope at none:2
[6] eval at .\boot.jl:331 [inlined]
[7] eval(::Expr) at .\client.jl:467
[8] top-level scope at .\none:3
in expression starting at C:\Users\Sannan.julia\packages\Yeppp\LvmJH\src\Yeppp.jl:8
ERROR: LoadError: Failed to precompile Yeppp [6310b701-1812-5374-a82f-9f6f2d54a40a] to
C:\Users\Sannan.julia\compiled\v1.5\Yeppp\xtx1U_lKfcn.ji.
Stacktrace:
[1] error(::String) at .\error.jl:33
[2] compilecache(::Base.PkgId, ::String) at .\loading.jl:1290
[3] _require(::Base.PkgId) at .\loading.jl:1030
[4] require(::Base.PkgId) at .\loading.jl:928
[5] require(::Module, ::Symbol) at .\loading.jl:923
[6] include(::Function, ::Module, ::String) at .\Base.jl:380
[7] include(::Module, ::String) at .\Base.jl:368
[8] top-level scope at none:2
[9] eval at .\boot.jl:331 [inlined]
[10] eval(::Expr) at .\client.jl:467
[11] top-level scope at .\none:3
in expression starting at C:\Users\Sannan.julia\packages\LowLevelParticleFilters\nvB4E\src\LowLevelParticleFilters.jl:7
ERROR: Failed to precompile LowLevelParticleFilters [d9d29d28-c116-5dba-9239-57a5fe23875b] to C:\Users\Sannan.julia\compiled\v1.5\LowLevelParticleFilters\m71p8_JHcjR.ji.
Stacktrace:
[1] error(::String) at .\error.jl:33
[2] compilecache(::Base.PkgId, ::String) at .\loading.jl:1290
[3] _require(::Base.PkgId) at .\loading.jl:1030
[4] require(::Base.PkgId) at .\loading.jl:928
[5] require(::Module, ::Symbol) at .\loading.jl:923

YEPPP ERROR 404 NOT FOUND

julia> Pkg.build("Yeppp")
Building Yeppp → C:\Users\Sannan\.julia\packages\Yeppp\LvmJH\deps\build.log
┌ Error: Error building Yeppp:
│ Exception calling "DownloadFile" with "2" argument(s): "The remote server returned an
│ error: (404) Not Found."
│ At line:1 char:82
│ + ... pe]::Tls12; (new-object net.webclient).DownloadFile("http://bitbucket ...
│ + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
│ + CategoryInfo : NotSpecified: (:) [], MethodInvocationException
│ + FullyQualifiedErrorId : WebException
│
│ [ Info: Attempting to create directory C:\Users\Sannan.julia\packages\Yeppp\LvmJH\deps\downloads
│ [ Info: Directory C:\Users\Sannan.julia\packages\Yeppp\LvmJH\deps\downloads already
exists
│ [ Info: Downloading file http://bitbucket.org/MDukhan/yeppp/downloads/yeppp-1.0.0.tar.bz2
│ ERROR: LoadError: failed process: Process('C:\Windows\System32\WindowsPowerShell\v1.0\powershell' -NoProfile -Command '[Net.ServicePointManager]::SecurityProtocol = [Net.SecurityProtocolType]::Tls12; (new-object net.webclient).DownloadFile("http://bitbucket.org/MDukhan/yeppp/downloads/yeppp-1.0.0.tar.bz2", "C:\Users\Sannan\.julia\packages\Yeppp\LvmJH\deps\downloads\yeppp-1.0.0.tar.bz2")', ProcessExited(1)) [1]
│
│ Stacktrace:
│ [1] pipeline_error at .\process.jl:525 [inlined]
│ [2] run(::Cmd; wait::Bool) at .\process.jl:440
│ [3] run(::Cmd) at .\process.jl:438
│ [4] run(::BinDeps.SynchronousStepCollection) at C:\Users\Sannan.julia\packages\BinDeps\ZEval\src\BinDeps.jl:521
│ [5] run(::FileRule) at C:\Users\Sannan.julia\packages\BinDeps\ZEval\src\BinDeps.jl:483
│ [6] run(::BinDeps.SynchronousStepCollection) at C:\Users\Sannan.julia\packages\BinDeps\ZEval\src\BinDeps.jl:521
│ [7] satisfy!(::BinDeps.LibraryDependency, ::Array{DataType,1}) at C:\Users\Sannan.julia\packages\BinDeps\ZEval\src\dependencies.jl:944
│ [8] satisfy!(::BinDeps.LibraryDependency) at C:\Users\Sannan.julia\packages\BinDeps\ZEval\src\dependencies.jl:922
│ [9] top-level scope at C:\Users\Sannan.julia\packages\BinDeps\ZEval\src\dependencies.jl:977
│ [10] include(::String) at .\client.jl:457
│ [11] top-level scope at none:5
│ in expression starting at C:\Users\Sannan.julia\packages\Yeppp\LvmJH\deps\build.jl:14
└ @ Pkg.Operations C:\Users\Sannan\AppData\Local\Programs\Julia 1.5.0\share\julia\stdlib\v1.5\Pkg\src\Operations.jl:942

Thank you!

Typo in DAEUnscentedKalmanFilter docstring

At the following line:

LowLevelParticleFilters.jl/src/ukf.jl

Line 190 in 3686062

 the measurements may be functions of both differential states `x` and algebraic variables `z`. Note, the actual dynamcis and measurement functions stored in the internal `ukf` should have signatures `(xz, u, p, t)`, i.e., they take the combined state containing both `x` and `z` in a single vector as dictated by the function `build_xz`. It is only the function `g` that is assumed to actually have the signature `g(x,z,u,p,t)`. 

"Note, the actual dynamcis" should be "Note, the actual dynamics".

If you think, I can open a PR to fix that.

Thanks

Struct of Arrays and BLAS2->3 optimization

Many operations can be sped up by considering batch evaluation of f/g/dg, i.e.,

A*x[i] + B*u[i] -> A*X + B*U

See LTI implementation here

Maybe an option that calls dynamics with the full particle vector would make this possible.

Vector{SVector} has the same layout as X/U above, just reinterpret.
Figure out how to handle the combined resampling x[i] = f(xp[j[i]], u, t, noise). Maybe by custom copyto!.

Such reinterpreted arrays can be modified inplace

julia> a = [SVector(1.0, 2.0) for _ in 1:3]
3-element Array{SArray{Tuple{2},Float64,1,2},1}:
 [1.0, 2.0]
 [1.0, 2.0]
 [1.0, 2.0]

julia> b = reinterpret(Float64, a)
6-element reinterpret(Float64, ::Array{SArray{Tuple{2},Float64,1,2},1}):
 1.0
 2.0
 1.0
 2.0
 1.0
 2.0

julia> b .+= 3
6-element reinterpret(Float64, ::Array{SArray{Tuple{2},Float64,1,2},1}):
 4.0
 5.0
 4.0
 5.0
 4.0
 5.0

LowLevelParticleFilters won't compile on Julia 1.8/macOS

Using julia 1.8.5 installed via juliaup on macOS 13.2.1

    julia> Pkg.precompile()
    Precompiling project...
      ✗ LowLevelParticleFilters
      ✗ ControlSystemIdentification
      0 dependencies successfully precompiled in 9 seconds. 305 already precompiled.
    
    ERROR: The following 2 direct dependencies failed to precompile:
    
    LowLevelParticleFilters [d9d29d28-c116-5dba-9239-57a5fe23875b]
    
    Failed to precompile LowLevelParticleFilters [d9d29d28-c116-5dba-9239-57a5fe23875b] to /Users/glenhenshaw/.julia/compiled/v1.8/LowLevelParticleFilters/jl_TVYZZ2.
    ERROR: LoadError: UndefVarError: SymbolicIndexingInterface not defined
    Stacktrace:
     [1] include(mod::Module, _path::String)
       @ Base ./Base.jl:419
     [2] include(x::String)
       @ LowLevelParticleFilters ~/.julia/packages/LowLevelParticleFilters/neZRi/src/LowLevelParticleFilters.jl:1
     [3] top-level scope
       @ ~/.julia/packages/LowLevelParticleFilters/neZRi/src/LowLevelParticleFilters.jl:28
     [4] include
       @ ./Base.jl:419 [inlined]
     [5] include_package_for_output(pkg::Base.PkgId, input::String, depot_path::Vector{String}, dl_load_path::Vector{String}, load_path::Vector{String}, concrete_deps::Vector{Pair{Base.PkgId, UInt64}}, source::Nothing)
       @ Base ./loading.jl:1554
     [6] top-level scope
       @ stdin:1
    in expression starting at /Users/glenhenshaw/.julia/packages/LowLevelParticleFilters/neZRi/src/PFtypes.jl:4
    in expression starting at /Users/glenhenshaw/.julia/packages/LowLevelParticleFilters/neZRi/src/LowLevelParticleFilters.jl:1
    in expression starting at stdin:1
    
    ControlSystemIdentification [3abffc1c-5106-53b7-b354-a47bfc086282]

Get rid of xprev

Almost never used but causes a lot of allocations. Get rid of permuting propagate_particles and have inplace resample instead.

Support for non-uniform observations?

Hello!

Before I try to call an external model in C++, I wanted to know if there is a support for non-uniform observations? For example, here's our constraint with a hydrological model. We want to assimilate snow data that are gathered "randomly" by human.

Our observations are not "periodic" : most years, but not always, observations are available on February 15th, March 1st-15th, April 1st
Rest of the year, there are no observations
Some year, we have continuous obs between February and May

Cheers|

Store exp(weights)

exp(w) is one of the most expensive computations, and it is done several times in

logsumexp
shouldresample
resample

Incremental precomiplation fatally broken

After updating all my packages, I got a warning from this package:

 WARNING: Method definition parameters(Any) in module SymbolicIndexingInterface at C:\Users\user\.julia\packages\SymbolicIndexingInterface\OL4Nh\src\interface.jl:71 overwritten in module LowLevelParticleFilters at C:\Users\user\.julia\packages\LowLevelParticleFilters\lSWeO\src\PFtypes.jl:9.
│    ** incremental compilation may be fatally broken for this module **

add_noise! assumes additive noise

Either force user to implement add_noise and save it in the APF struct or do a second round of propagate_particles with noise (expensive but general).

Register

@JuliaRegistrator register()

Installing this package is downgrading many other packages

julia> versioninfo()
Julia Version 1.10.0
Commit 3120989f39b (2023-12-25 18:01 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 32 × AMD Ryzen 9 7950X 16-Core Processor
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, znver3)
  Threads: 1 on 32 virtual cores
Environment:
  LD_LIBRARY_PATH = /lib:/usr/lib:/usr/local/lib

MWE:

mkdir test
cd test
julia --project="."

and then in Julia:

]
add ModelingToolkit

This prints two errors, you can ignore them. Restart Julia, and then:

]
st
add LowLevelParticleFilters

gives the output:

(test) pkg> st
Status `~/repos/test/Project.toml`
  [961ee093] ModelingToolkit v8.75.0

(test) pkg> add LowLevelParticleFilters
   Resolving package versions...
   Installed SeeToDee ──────────────── v1.2.0
   Installed FastGaussQuadrature ───── v0.5.1
   Installed LowLevelParticleFilters ─ v3.6.1
   Installed ModelingToolkit ───────── v8.73.2
   Installed OrdinaryDiffEq ────────── v6.59.3
    Updating `~/repos/test/Project.toml`
  [d9d29d28] + LowLevelParticleFilters v3.6.1
⌃ [961ee093] ↓ ModelingToolkit v8.75.0 ⇒ v8.73.2
    Updating `~/repos/test/Manifest.toml`
⌅ [c3fe647b] + AbstractAlgebra v0.34.7
  [7d9f7c33] + Accessors v0.1.35
  [4c555306] - ArrayLayouts v1.5.2
  [a33af91c] + CompositionsBase v0.1.2
  [8bb1440f] + DelimitedFiles v1.9.1
⌃ [2b5f629d] ↓ DiffEqBase v6.146.0 ⇒ v6.145.2
⌃ [459566f4] ↓ DiffEqCallbacks v2.36.1 ⇒ v2.35.0
⌅ [442a2c76] + FastGaussQuadrature v0.5.1
⌅ [0b43b601] + Groebner v0.5.1
⌅ [d5909c97] + GroupsCore v0.4.2
  [18e54dd8] + IntegerMathUtils v0.1.2
  [3587e190] + InverseFunctions v0.1.12
  [73f95e8e] + LatticeRules v0.0.1
  [5078a376] - LazyArrays v1.8.3
  [d9d29d28] + LowLevelParticleFilters v3.6.1
  [a3b82374] - MatrixFactorizations v2.1.0
  [bb5d69b7] - MaybeInplace v0.1.1
⌃ [961ee093] ↓ ModelingToolkit v8.75.0 ⇒ v8.73.2
⌅ [8913a72c] ↓ NonlinearSolve v3.4.0 ⇒ v2.8.2
⌃ [1dea7af3] ↓ OrdinaryDiffEq v6.69.0 ⇒ v6.59.3
  [08abe8d2] + PrettyTables v2.3.1
  [27ebfcd6] + Primes v0.5.5
  [8a4e6c94] + QuasiMonteCarlo v0.3.3
  [fb686558] + RandomExtensions v0.4.4
⌅ [731186ca] ↓ RecursiveArrayTools v3.5.4 ⇒ v2.38.10
  [fdea26ae] + SIMD v3.4.6
⌃ [0bca4576] ↓ SciMLBase v2.20.0 ⇒ v2.10.0
  [e9a6253c] + SciMLNLSolve v0.1.9
⌃ [1c904df7] + SeeToDee v1.2.0
⌅ [727e6d20] ↓ SimpleNonlinearSolve v1.3.1 ⇒ v0.1.25
  [ed01d8cd] + Sobol v1.5.0
  [892a3eda] + StringManipulation v0.3.4
⌅ [2efcf032] ↓ SymbolicIndexingInterface v0.3.3 ⇒ v0.2.2
⌃ [d1185830] ↓ SymbolicUtils v1.5.0 ⇒ v1.4.0
⌃ [0c5d862f] ↓ Symbolics v5.16.1 ⇒ v5.11.0
        Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated -m`
Precompiling project...
  32 dependencies successfully precompiled in 111 seconds. 193 already precompiled.

(test) pkg>

Expected result: Installing ModelingToolkit and then LowLevelParticleFilters is NOT downgrading any packages.

Does `correct!()` support `missing` samples as explained in docs?

Hi, this issue pertains to a prior issue and this section of the docs

I was digging through the code to understand whether correct!() could be modified to bypass updating the weights of particles in the case that the observation y[i] == missing, i.e. regular sampling rate but with sparse dropout of measurements. Based on the prior issue I was led to believe that this it is not the case that correct!() knows about missing data and one would have to write a manual loop to handle this case (as the documentation suggests).

However, the correct!() function calls measurement_equation!() which has a short circuit evaluation any(ismissing.(y)) && return w

To me, this looks like it would mean that functions such as correct!() does know to "ignore" missing observations, and therefore higher-level loops like forward_trajectorry() indeed handle missing data already. If how I am understanding it is truly what's going on, then could the documentation be updated to reflect this (i.e. the loops used by forward_trajectorry() and the like do not apply corrections if the observation for that time step is missing)? Else, I'd like to know whether the interpretation of the short circuit should be.

Thanks!

Compatibility to ModelingToolkit.jl?

Hello Fredrik!

Thanks for the great package! :)

I want to icorporate the package for nonlinear state estimation. I have built an ODE model using MTK. Is there a straightforward way to use the estimators together with MTK models (ODESystem type)?

I assume I could try to convert my ODESystem to a ODEProblem from DifferentialEquations.jl and obtain the state transition function by the f-field. Or is there another way to automatically "de-modelingtoolkitize" the MTK model from symbolic to a numeric function in order to use it with the estimators?

Otherwise I would have to rewrite the model to a numeric function, right?

What is your thought on that?

loglik bug for ParticleFilter with resampling

The standard ParticleFilter returns incorrect loglikelihood unless resampling is carried out every step. This indicates that there is some issue with how maxw is propagated since this is reset to 0 upon resampling.

v3.0

parameter argument
Store parameters in filters
DAE estimation (Nonlinear State Estimation of Differential Algebraic Systems, Mandela et al.)
FilteringProblem interface, (smoothing etc. as well)
Time varying covariance matrices
kalman matrices can be functions of parameters and time (state and control as well)
fixed-gain observer
remove some methods that were added for compat when adding u to correct!
result structures

Next major version

Time instead of index (or in addition to?)
initial state distribution part of filter or argument to forward_trajectory etc.?
Filters to store Ts. If default Ts::Int = 1, we can have t be time without breaking the interface.

Remove StatsPlots dependency

StatsPlots is used only for diagnostics, it's not necessarily needed for the package to function correctly. Combined with its relatively large size, I suggest to remove it and only keep it in the examples files with a note about it being needed for the examples to function.

Improve documentation

There is an "advanced tutorial" section in the manual, but I miss a "basic tutorial".

I could add the example we discussed here https://discourse.julialang.org/t/discrete-kalman-filter/109014 as pull request to the documentation, also extend it a bit as an example for amending a dynamical system to estimate an additional state (Ta)...

What do you think? And if you think this is a good idea, where in the documentation should I add it?

Passing dynamics noise density makes stochastic dynamics less expressive

I really like the level of flexibility this package affords compared to other particle filtering packages in Julia, so first of all, thank you for developing and maintaining it. I had a quick comment about the source code, where the propagation of the dynamics model occurs. The dynamics model has an option for injecting noise into the dynamics (set to false as default), so shouldn't this be utilized in the propagation step? The following code makes it seem like you can only make the dynamic noise be additive:

Base.@propagate_inbounds function propagate_particles!(pf::ParticleFilter,u,j::Vector{Int}, p, t::Int, d::Union{Nothing, Distributions.Sampleable}=pf.dynamics_density)
    f = dynamics(pf)
    s = state(pf)
    x,xp = s.x, s.xprev
    VecT = eltype(s.x)
    D = length(VecT)
    noise = zeros(D)
    if d === nothing
        for i = eachindex(x)
            x[i] =  f(xp[j[i]], u, p, t)
        end
    else
        for i = eachindex(x)
            x[i] =  f(xp[j[i]], u, p, t) + VecT(rand!(pf.rng, d, noise))
        end
    end
    x
end

I was thinking something like this would make more sense to let the dynamics be arbitrarily stochastic:

...
    if d === nothing
        for i = eachindex(x)
            x[i] =  f(xp[j[i]], u, p, t)
        end
    else
        for i = eachindex(x)
            x[i] =  f(xp[j[i]], u, p, t, true)
        end
    end
...

To be quite frank, I would think you don't even need to pass the dynamics noise density into the filter to initialize it since it'll be redundant (just like we don't pass the observation noise density), and we need only pass the initial condition density.