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

**URL:** <https://discourse.julialang.org/t/rmse-vague-in-2d-particle-filtering-with-lowlevelparticlefilters-jl/131042>\
**Category:** Modelling & Simulations\
**Tags:** question\
**Created:** [July 25, 2025, 5:06am UTC](https://discourse.julialang.org/t/rmse-vague-in-2d-particle-filtering-with-lowlevelparticlefilters-jl/131042 "2025-07-25T05:06:33Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![Tanya123](https://avatars.discourse-cdn.com/v4/letter/t/a8b319/32.png) [@Tanya123](https://discourse.julialang.org/u/Tanya123)\
**Post date:** [July 25, 2025, 5:06am UTC](https://discourse.julialang.org/t/rmse-vague-in-2d-particle-filtering-with-lowlevelparticlefilters-jl/131042/1 "2025-07-25T05:06:33Z")

</div>

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 :

```julia-auto
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](https://global.discourse-cdn.com/julialang/original/3X/7/e/7efb2401dbbd31835489a8308f878f585a8ccac3.png)

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.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [July 25, 2025, 6:37am UTC](https://discourse.julialang.org/t/rmse-vague-in-2d-particle-filtering-with-lowlevelparticlefilters-jl/131042/2 "2025-07-25T06:37:13Z")

</div>

See answer here

> [@Particle Filtering methods](https://discourse.julialang.org/t/particle-filtering-methods/130828/10):
>
> It depends on how large the uncertainty is in the initial distribution compared to the stationary distribution. If you push that example much further, it might converge. You can compute the stationary covariance matrix by solving a Riccati equation with the A jacobian, see [Influence of sample rate on performance · LowLevelParticleFilters Documentation](https://baggepinnen.github.io/LowLevelParticleFilters.jl/stable/sample_rate/) for some details or “Asymptotic form” in [Kalman filter - Wikipedia](https://en.wikipedia.org/wiki/Kalman_filter#Asymptotic_form) using MatrixEquations A = Ajac(...) C = Cjac(...) R\_stationary, \_ = ared(A, C'…

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