diff --git a/docs/src/beetle_example.md b/docs/src/beetle_example.md index 9fa5708f..d29bdcce 100644 --- a/docs/src/beetle_example.md +++ b/docs/src/beetle_example.md @@ -108,8 +108,8 @@ plot!(xyt[:,1],xyt[:,2], c=:red, lab="measurement") ``` as well as the angle state variable (we subsample the particles to not get sluggish plots) ```@example beetle -fig2 = scatter(to1series(ϕ.(x)'[:,1:5:end])..., m=(:black, 0.03, 2), lab="", size=(500,300), format=:png) -plot!(identity.(xh[:,4]), lab="Filtered angle", legend=:topleft, ylims=(-30, 70)) +fig2 = scatter(to1series(ϕ.(x)'[:,1:5:end])..., m=(:black, 0.03, 2), lab="", size=(500,300), format=:png, dpi=100) +plot!(identity.(xh[:,4]), lab="Filtered angle", legend=:topleft, ylims=(-30, 70), format=:png) ``` The particle plot above indicate that the posterior is multimodal. This phenomenon arises due to the simple model that uses an angle that is allowed to leave the interval ``0-2\pi` rad. In this example, we are not interested in the angle, but rather when the beetle switches mode. The filtering distribution above gives a hint at when this happens, but we will not plot the mode trajectory until we have explored smoothing as well. diff --git a/docs/src/index.md b/docs/src/index.md index 59d8624d..ca1840e6 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -289,7 +289,7 @@ nothing # hide We now create the [`AdvancedParticleFilter`](@ref) and use it in the same way as the other filters: ```@example lingauss -apf = AdvancedParticleFilter(N, dynamics, measurement, measurement_likelihood, df, d0) +apf = AdvancedParticleFilter(N, dynamics, measurement, measurement_likelihood, d0) sol = forward_trajectory(apf, u, y, p) ``` diff --git a/src/PFtypes.jl b/src/PFtypes.jl index 1d8ddf5b..21559d27 100644 --- a/src/PFtypes.jl +++ b/src/PFtypes.jl @@ -164,7 +164,7 @@ end dynamics::FT measurement::GT measurement_likelihood::GLT - dynamics_density::FDT = Normal() + dynamics_density::FDT = MissingDistribution() initial_density::IDT resample_threshold::Float64 = 0.5 resampling_strategy::RST = ResampleSystematic @@ -175,7 +175,7 @@ end """ - AdvancedParticleFilter(N::Integer, dynamics::Function, measurement::Function, measurement_likelihood, dynamics_density, initial_density; p = SciMLBase.NullParameters(), threads = false, kwargs...) + AdvancedParticleFilter(N::Integer, dynamics::Function, measurement::Function, measurement_likelihood, initial_density; p = SciMLBase.NullParameters(), threads = false, kwargs...) This type represents a standard particle filter but affords extra flexibility compared to the [`ParticleFilter`](@ref) type, e.g., non-additive noise in the dynamics and measurement functions. @@ -186,7 +186,6 @@ See the docs for more information: https://baggepinnen.github.io/LowLevelParticl - `dynamics`: A discrete-time dynamics function `(x, u, p, t, noise=false) -> x⁺`. It's important that the `noise` argument defaults to `false`. - `measurement`: A measurement function `(x, u, p, t, noise=false) -> y`. It's important that the `noise` argument defaults to `false`. - `measurement_likelihood`: A function `(x, u, y, p, t)->logl` to evaluate the log-likelihood of a measurement. -- `dynamics_density`: This field is not used by the advanced filter and can be set to `nothing`. - `initial_density`: The distribution of the initial state. - `threads`: use threads to propagate particles in parallel. Only activate this if your dynamics is thread-safe. `SeeToDee.SimpleColloc` is not thread-safe by default due to the use of internal caches, but `SeeToDee.Rk4` is. """ @@ -202,6 +201,9 @@ function AdvancedParticleFilter(N::Integer, dynamics::Function, measurement::Fun initial_density, kwargs...) end +# Constructor without dynamics_density which is not used except for in smoothing +AdvancedParticleFilter(N::Integer, dynamics::Function, measurement::Function, measurement_likelihood, initial_density; kwargs...) = AdvancedParticleFilter(N, dynamics, measurement, measurement_likelihood, MissingDistribution(), initial_density; kwargs...) + Base.@propagate_inbounds function measurement_equation!(pf::AbstractParticleFilter, u, y, p, t, w = weights(pf)) g = measurement_likelihood(pf) diff --git a/src/smoothing.jl b/src/smoothing.jl index 9bae2d8d..71a301c6 100644 --- a/src/smoothing.jl +++ b/src/smoothing.jl @@ -43,6 +43,7 @@ function smooth(pf::AbstractParticleFilter, xf, wf, wef, ll, M, u, y, p=paramete N = num_particles(pf) f = dynamics(pf) df = dynamics_density(pf) + df isa MissingDistribution && error("In order to perform smoothing, you must provide a noise density for the dynamics. Use the constructor `AdvancedParticleFilter(N, dynamics, measurement, measurement_likelihood, dynamics_density, initial_density; kwargs...)` to do this") @assert M <= N "Must extend cache size of bins and j to allow this" xb = Array{particletype(pf)}(undef,M,T) j = resample(ResampleSystematic, wef[:,T], M) diff --git a/src/utils.jl b/src/utils.jl index a81c9cea..2b7ab2ed 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -115,6 +115,8 @@ Mixed value support indicates that the distribution is a mix of continuous and d """ struct Mixed <: ValueSupport end +struct MissingDistribution <: MultivariateDistribution{Mixed} end + """ TupleProduct(v::NTuple{N,UnivariateDistribution}) diff --git a/test/runtests.jl b/test/runtests.jl index 656a0d51..395e4e4f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -292,7 +292,7 @@ end T = 200 # Number of time steps N = 500 Random.seed!(0) - apf = AdvancedParticleFilter(N, dynamics, measurement, measurement_likelihood, df, d0, threads=false) + apf = AdvancedParticleFilter(N, dynamics, measurement, measurement_likelihood, d0, threads=false) x,u,y = LowLevelParticleFilters.simulate(apf,T,du) # Simuate trajectory using the model in the filter x,u,y = tosvec.((x,u,y)) @time resapf,ll = mean_trajectory(apf, u, y) @@ -302,7 +302,7 @@ end # With threads Random.seed!(0) - apf = AdvancedParticleFilter(N, dynamics, measurement, measurement_likelihood, df, d0, threads=true) + apf = AdvancedParticleFilter(N, dynamics, measurement, measurement_likelihood, d0, threads=true) x,u,y = LowLevelParticleFilters.simulate(apf,T,du) # Simuate trajectory using the model in the filter x,u,y = tosvec.((x,u,y)) @time resapf,ll = mean_trajectory(apf, u, y)