RMSE vague in 2D particle filtering with LowLevelParticleFilters.jl?

Specific Domains Modelling & Simulations
1

Hi all,

I’m applying particle filtering to a 2D linear advection problem (60x60 grid, 3600 states) using LowLevelParticleFilters.jl with:

*Particles: 10,000
*Process noise: 0.01
*Measurement noise: 0.1 at 100 random sensor points
*Dynamics: linear advection with periodic BCs

Code :

using Pkg
Pkg.add("LowLevelParticleFilters")
Pkg.add("LinearAlgebra")
Pkg.add("Random")
Pkg.add("Statistics")
Pkg.add("Distributions")
Pkg.add("Plots")


using LowLevelParticleFilters
using LinearAlgebra
using Random
using Statistics
using Distributions
using Plots

function main()
    # Grid parameters
    Nx, Ny = 60, 60
    N = Nx * Ny
    T = 100            # number of time steps
    Np = 10000           # number of particles

    # Wavenumber
    w = 2π / 1.0

    # Initial field (sinusoidal bump)
    function initial_field()
        x = LinRange(0, 1.0, Nx)
        y = LinRange(0, 1.0, Ny)
        [1000 + 1000*sin(w*xi) * sin(w*yi) for xi in x, yi in y] |> vec
    end

    # Dynamics: linear advection to the right
    function dynamics(x, u, p, t)
        β = reshape(x, Nx, Ny)
        u_x = 1.0
        dx = 2π / Nx
        dt = 0.2 * dx / u_x

        β_left = circshift(β, (0, -1))
        dbdx = (β - β_left) / dx
        β_new = β .- u_x * dt .* dbdx
        return vec(β_new)
    end

    # Measurement: sample at random locations
    sensor_indices = rand(1:N, 100)
    function measurement(x, u, p, t)
        x[sensor_indices]
    end

    # Noise models
    process_noise = MvNormal(zeros(N), 0.01^2 * I(N))
    measurement_noise = MvNormal(zeros(length(sensor_indices)), 0.1^2 * I(length(sensor_indices)))
    initial_state_dist = MvNormal(initial_field(), 0.1^2 * I(N))

    # Initialize PF
    pf = ParticleFilter(Np, dynamics, measurement, process_noise, measurement_noise, initial_state_dist)

    # Generate true trajectory & observations
    x_true = initial_field()
    true_states = Vector{Vector{Float64}}(undef, T)
    observations = Vector{Vector{Float64}}(undef, T)

    for t in 1:T
        x_true = dynamics(x_true, nothing, nothing, t) + rand(process_noise)
        y = measurement(x_true, nothing, nothing, t) + rand(measurement_noise)
        true_states[t] = x_true
        observations[t] = y
    end

    # Run PF
    estimates = Vector{Vector{Float64}}(undef, T)
    for t in 1:T
        @show t
        pf(nothing, observations[t])
        estimates[t] = weighted_mean(pf)
    end

    # RMSE over time
    rmse_history = zeros(T)
    for t in 1:T
        err = estimates[t] - true_states[t]
        rmse_history[t] = sqrt(mean(err.^2))
    end

    # RMSE Plot
    plt_rmse = plot(
        1:T, rmse_history,
        xlabel = "Time step",
        ylabel = "RMSE",
        title = "RMSE Evolution over Time",
        linewidth = 2,
        marker = :circle,
        legend = false
    )
    savefig(plt_rmse, "rmse_particle_filter.png")
    display(plt_rmse)

    println("✅ Particle Filter completed and RMSE plot saved as rmse_particle_filter.png")
end

main()

(Thanks to @baggepinnen for guidance on the dynamics and for sharing EKF code that I modified for particle filtering in this implementation.)

The RMSE over time is vague and does not show clear convergence, and I am unable to reason about it.

Image:

image
image600×400 23.8 KB

Is this implementation fundamentally correct, or is there something structurally wrong in the code? I understand the limitations of particle filtering in very high-dimensional settings, but I would like to confirm whether the code is functioning as intended (even if performance is poor), or if there are critical mistakes I should address first.

Any help or comment on this will be really appreciated.

2

See answer here

particle filters add additional monte-carlo approximation error that makes the analysis slightly harder.