How many particles, and how wrong can your first guess be?

Three posts ago I derived the Kalman filter from a voltmeter. Two posts ago its extended cousin locked onto the wrong range and stayed there. Last post a cloud of four thousand weighted guesses held two stories at once (close-and-slow, far-and-fast) for exactly as long as the lighthouse’s one bearing gave it no reason to pick. That is the whole arc: 1-D → nonlinear → multi-modal, each post breaking the previous one’s assumption. This one closes it by asking the two questions nobody answers with a derivation: how many particles is enough, and how wrong can the first guess be before it matters. §5.5–5.7 and §6.5–6.7 of the report (Experiments 5, 6 and 7) answer both, with a third experiment thrown in that answers a question I didn’t expect: what happens to the multi-modality itself when you add a second kind of measurement.

The code for this post is Experement_5.py (and its loader), Experement_6.py (and its loader), Experement_7.py, Problems.py::TestProblem4 and two more functions in Utils.py. The widget is new work, not a port of an experiment script: a Rust crate, pf-wasm, compiled to WebAssembly and run in a Web Worker, because the honest budget for “run the sweep live” is tens of millions of particle-updates: the one place in this four-post arc where that toolchain earns its keep rather than being a distraction.

Question one: how many particles?

§5.6 sweeps $N_s$ : the report says “from $N_s=100$ to $N_s=10000$ using steps of 100”; the script that actually produced the figure, Experement_6.py line 45, loops range(200, 10200, 200), 50 points from 200 to 10000. Small discrepancy, but worth noting since I’m about to lean on the code, not the prose, for the widget’s own point count. Everything else is fixed at dataset B’s usual values (the same $\mathbf{Q}$ , $\mathbf{R}$ , $K=0.2$ as post 19) except the initial guess, which here is the near cluster, not the true start:

initP = np.array([[0.5**2,       0,      0,       0],
                  [0     , 0.005**2,      0,       0],
                  [0     ,        0, 0.3**2,       0],
                  [0     ,        0,      0, 0.01**2]])
initX = [0.0, 0.0, 0.4, -0.05]
...
runs = 20
def evaluateAllFiltersPartial(numberOfParticals, currX):
    return evaluatePFFilterLessMemory(problem, numSamples, X0Real, numberOfParticals, K,
                               R, Q, initP, initX, [0,1,2,3])

Experement_6.py lines 11–33 (trimmed). evaluatePFFilterLessMemory (Utils.py lines 100–114) only ever builds one filter, ParticleFilterSIRWithJitter (SIR, always resampling, jitter on) and runs it 20 times per $N_s$ . Fig. 12, page 22:

The report’s own advice, §7 p.24, is the whole reason this post has a live sweep instead of a screenshot of one: picking $N_s$ by “however many the CPU can afford” wastes resources, and picking it from a general rule of thumb gets the wrong answer for your $f$ , $h$ and dimensionality. The right move is to run the algorithm a few times at increasing $N_s$ , plot error against $N_s$ , and read the knee off the curve. So: press play below rather than read mine off a report page.

Two stacked log-scale charts, EKF mean squared error above and particle-filter mean squared error below, plotted against the sweep axis with a shaded interquartile band around a median line. A dashed marker on the lower chart flags the smallest particle count within 15 percent of the best error seen. — InteractiveRun the sweep yourself

I ran it myself at $N_s$ from 100 to 10000 (GPF with jitter, 20 runs a point, matching Experement_6.py’s own runs = 20): the median MSE drops from about $1.0\times10^{-2}$ at $N_s=100$ to $4$ – $7\times10^{-3}$ by $N_s\approx1000$ and stays in that band out to 10000, the same decaying-then-flat shape as Fig. 12, noisier because 20 runs a point is what the report itself budgeted, not because the widget cut corners. The knee marker (dashed, “smallest $N_s$ within 15% of the best seen”) lands around 1000–2000 in my runs, a little earlier than the report’s “sufficient by 4000”, both are reading noise off the same shape, and that is rather the point: your curve, for your problem, is the one to trust, and this one visibly moves as more points stream in.

Question two: how wrong can the first guess be?

§5.5 asks the complementary question: fix $N_s$ , and instead push the initial estimate away from the truth. Experement_5.py lines 20–46:

X0Real = np.array([-0.05, 0.001, 0.7, -0.055])
weightings = [-15.0, 8.0, 15.0, -8.0]
...
def evaluateAllFiltersPartial(initX, currX):
    initP = np.eye(4)*(1e-6+(X0Real-initX)**2)
    return evaluateImportantAllFiltersLessMemory(problem, numSamples, X0Real, numberOfParticals, K,
                               R, Q, initP, initX, [0,1,2,3])
...
for axis in range(4):
    for percentageOfActual in np.linspace(0,1,100):
        X0RealValue = X0Real.copy()
        X0RealValue[axis] += weightings[axis]*percentageOfActual

Read that carefully, because it is a nice piece of engineering: the true start X0Real never moves: only the guess X0RealValue walks away from it, one state axis at a time, up to $\pm15$ in position or $\pm8$ in velocity. And the guess’s own uncertainty, initP, grows to match: $\mathbf{P}_0 = 10^{-6}\mathbf{I} + (\mathbf{x}^{\text{real}}_0-\mathbf{x}^{kf}_0)^2$ . The filter is told, correctly, “I might be off by this much”, which for a particle cloud means the true state sits about one standard deviation from the sampled mean at every point on the sweep, always recoverable in principle. Whether it is recovered in practice is the experiment.

Figs. 10 and 11, page 21, split the two filters onto separate charts, because their axes don’t share a scale: the EKF’s climbs to $10^4$ – $10^8$ while the PF’s stays at $10^{-4}$ – $10^{-2}$ , four to ten orders of magnitude apart. My own re-run (widget defaults, GPF with jitter, $N_s=2000$ , 50 runs a point: the plan’s honest budget for a browser sweep, half the report’s $N_s=4000$ and a twentieth of its 1000 runs) lands on the same shape and very nearly the same scale:

Axis (offset at full throw)	EKF MSE at offset 0	EKF MSE at full offset	PF MSE at full offset
$x_0$ ( $x$ , $\pm15$ )	$5.6\times10^{-4}$	$5.9\times10^4$	$6.6\times10^{-4}$
$x_1$ ( $\dot x$ , $\pm8$ )	$6.0\times10^{-4}$	$1.1\times10^6$	see below
$x_2$ ( $y$ , $\pm15$ )	$7.1\times10^{-4}$	$1.3\times10^5$	$4.7\times10^{-4}$
$x_3$ ( $\dot y$ , $\pm8$ )	$5.1\times10^{-4}$	$3.0\times10^8$	$3.3\times10^{-2}$

The particle filter earns its reputation on three of those four columns, barely moving while the EKF detonates. The fourth, $\dot x$ , is the honest exception, and it is worth showing rather than hiding.

The particle filter can get stuck too

At full throw on $x_1$ ( $\dot x$ ), my widget’s PF median MSE across 100 runs at $N_s=4000$ (the report’s own particle count) was $13.7$ , not the $10^{-3}$ Fig. 11’s axis suggests. Stepping through one run explains why: the true $x$ -velocity is $0.001$ , next to nothing, because this target’s motion is almost entirely along $y$ . Guess a $\dot x$ of $8$ instead and the cloud’s position estimate drifts by roughly 8 units a step: 24 steps is a drift of nearly 200, while the bearing sensor barely constrains $\dot x$ at all, since $x$ contributes little to $\operatorname{atan2}(y,x)$ while $y$ dominates. A handful of particles do start near the true velocity (it sits exactly one $\sigma$ from the sampled mean, by construction of initP), but if none of them survive the first resample, the whole cloud can lock onto a self-consistent, wrong trajectory, and once it has, the jitter that is supposed to rescue it (roughen, post 19’s $K\,E_j\,N_s^{-1/d}$ ) scales with the cloud’s current spread, which has just collapsed. A particle filter that has converged to the wrong answer looks exactly like one that has converged to the right one; nothing in the algorithm tells you which. This is the report’s own “particle filter can recover if a single particle survives in a low-probability region” (§7) read backwards: on this axis, over these seeds, often none did. It happens far less at the report’s full $N_s=4000$ /1000-run budget than it did in my faster, noisier one, but it is not zero, and a real engineer should see it happen at least once before trusting a particle filter blind. Press Run on the $\dot x$ ( $x_1$ ) axis a few times and watch for it.

Why: the report’s own explanation

§7’s closing paragraphs, rephrased: the further the initial guess sits from the truth, the faster the true object would have to be moving to produce the same sequence of bearings, so distance from the guess, along a velocity axis especially, doesn’t just add error, it manufactures alternative stories for the data to tell. That is literally the source of the bimodality post 19 built its whole argument on; Experiment 5 is the same phenomenon shown as a dial rather than a snapshot. And the flip side, worth remembering before the next section: an object that is close and moving fast has few alternative explanations: it’s already near the edge of what the sensor can distinguish, so the EKF does comparatively better the closer and faster the target is. Distance from truth is what buys the particle filter its advantage, not multi-modality as some intrinsic property of bearings-only tracking.

Question three (a bonus): what if you add a range measurement?

Experiment 7 changes the problem, not the initial guess. Dataset C (Problems.py::TestProblem4, lines 195–222) keeps the bearing but adds a second measurement, the squared range:

class TestProblem4(ProblemType):
    def __init__(self, r=0.005**2, q=0.001**2):
        self.n,self.m = 4,2
        ...
    def h(self, x,k):
        return np.array([np.nan_to_num(np.arctan2(x[2],x[0]))
        ,x[2]**2+x[0]**2])
    ...
    def Jh(self, x,k):
        x2py2=np.nan_to_num(1.0/(x[0]**2+x[2]**2))
        return np.array([[-x[2]*x2py2, 0, x[0]*x2py2, 0],
                         [2*x[0], 0, 2*x[2], 0]])

Squared range, not range, which is why the Jacobian’s second row is the clean $[2x_0,0,2x_2,0]$ rather than needing a square root anywhere. Experement_7.py sets $\mathbf{R}=\operatorname{diag}(0.005^2,0.0001)$ , keeps $N_s=4000$ , and runs all eight algorithms 1000 times. Table 3, page 22:

Algorithm	EKF	IEKF	SIS	SISwJ	GPF	GPFwJ	SIR	SIRwJ
Mean “RMSE” ( $\times10^{-3}$ )	5.43	0.31	0.89	1.38	0.14	0.23	0.16	0.44

The point isn’t the table, though: it’s what happened to the EKF compared with dataset B. Bearing-only, its mean MSE on this same starting point is $6.2$ (my re-run; the report’s own dataset-B baseline, §6.3, is $6.52$ , see post 19). Add the range measurement and it drops to $3.9\times10^{-3}$ , better than a thousand-fold. The particle filter improves too (mean MSE roughly $7\times10^{-3}\to7\times10^{-4}$ in my re-run, about $10\times$ ), but nowhere near as dramatically, because it didn’t need saving in the first place. The moral the plan for this post insists on, and I now believe: it is not that particle filters beat Kalman filters. It’s that a Gaussian can’t describe a multi-modal posterior, and a range measurement is one way to stop the posterior being multi-modal at all. Once it isn’t, the EKF is not just “usable again”: it goes back to being the better tool, because it is cheaper.

Why: cost

§7’s last real paragraph, rephrased, because it is the sentence I’d underline if I were reading this report cold: on a simple posterior, use the Kalman filter, because it is faster: its most expensive step is one matrix inversion (or a closed form, for small state dimensions), while a particle filter’s cost is $N_s$ times the per-particle work, run for every one of those particles. The saving grace is that per-particle work is embarrassingly parallel (a GPU or a many-core machine amortises it in a way a single matrix inversion cannot be split) and that particle filters spend their computation where the posterior actually has mass, which for a genuinely multi-modal problem, no fixed Gaussian ever can.

The crate: `pf-wasm`

The sweep above is real Monte Carlo, not a stored table: pick a point on an axis, and the browser draws 20–50 independent noisy datasets, runs an EKF and a full particle-filter cloud over 24 steps each, and repeats it across up to 40 points. At the top of that range (the offset sweep at $N_s=2000$ , 50 runs, 40 points) that’s on the order of $2000\times4\times24\times50\times40 \approx3.8\times10^8$ particle-state touches, not counting the per-particle weighing and resampling passes. Plain JS handles this fine at post 19’s scale (a few thousand particles at 60 fps); asked to do it 40 times over inside a progress bar, it is the difference between “a few seconds” and “the better part of a minute”, which is exactly the gap wasm/crates/pf-wasm/ exists to narrow.

The crate is self-contained (calib-wasm is the camera-calibration crate, gp-wasm the Gaussian-process one, nothing existing fit), follows example-wasm’s C-ABI convention (wasm_alloc/wasm_free verbatim, crate-type = ["cdylib"]), and has no dependencies, matching every other crate in this workspace: rng.rs (a bit-exact Rust port of the site’s world-sim seeded PRNG, sfc32 through splitmix32, verified against five next() and five normal() values recorded from the TypeScript class at seed 42), model.rs (dataset B/C, Problems.py::TestProblem3/TestProblem4), ekf.rs (ExtendedKalmanFilter.py, ported faithfully, including a bug I found while doing it, below), pf.rs (ParticalFilter.py’s SIS/GPF/SIR, log-domain weights, systematic resampling and the corrected jitter formula post 19 already established), and sweep.rs, which runs one point of the chart and reports mean/median/LQ/UQ for both filters. One exported function, sweep_point, does the whole inner loop per call so a Web Worker can call it once per sweep point and postMessage the result, which is how the chart streams in.

A bug in ExtendedKalmanFilter's own process noise

ExtendedKalmanFilter.__init__ (ExtendedKalmanFilter.py lines 21–28) takes r and q and stores self.R = r*np.eye(m), self.Q = q*np.eye(n). Every call site in Utils.py passes the actual matrices $\mathbf{R}$ and $\mathbf{Q}$ for r and q, not scalars. For $\mathbf{R}$ that’s harmless: $m=1$ (or 2), so R * eye(m) gives back R. But $\mathbf{Q}$ is the full correlated $4\times4$ process-noise matrix (the same one that generates the truth trajectory and the one the particle filter samples from) and Q * eye(4) is an elementwise product with the identity: it keeps the diagonal and zeroes every off-diagonal term. So every EKF result in this report, dataset A onward, was computed with $\mathbf{Q}_{\text{EKF}}=\operatorname{diag}(0.25q,q,0.25q,q)$ , the position–velocity correlation the report’s own §2 derives (the $\Delta t^3/2$ term of the white-noise-acceleration model) never actually reaches the filter’s covariance update. I didn’t fix it: pf-wasm’s ekf.rs reproduces q_ekf_diag on purpose, because the numbers above are supposed to be this filter’s numbers, not a hypothetically-correct one’s. I don’t know how large the effect is without re-running every earlier post’s EKF with the correlation restored, which I haven’t done, flagging it here rather than quantifying it.

cargo test -p pf-wasm has seven cases, chosen to catch exactly the classes of bug the sibling posts have found in this codebase: the PRNG matches the TypeScript class bit-for-bit at a known seed, not just itself across two instances; systematic resampling preserves the particle count and concentrates on a single heavy particle (I rig one particle to hold nearly all the weight and check it survives resampling at better than 90%); log-sum-exp agrees with the naive computation on well-scaled values and stays finite where the naive one underflows to zero (the whole reason scaleLogW exists, see post 19); and a full sweep point is bit-for-bit deterministic given a seed, run twice.

#[test]
fn log_sum_exp_matches_naive_when_naive_would_underflow() {
    let mut logw: Vec<f64> = vec![-800.0, -800.5, -799.2, -805.0];
    let naive: Vec<f64> = logw.iter().map(|v| v.exp()).collect();
    assert!(naive.iter().all(|&v| v == 0.0), "test setup: these should underflow");
    scale_log_w(&mut logw);
    let sum: f64 = logw.iter().map(|v| v.exp()).sum();
    assert!((sum - 1.0).abs() < 1e-9, "normalised weights should sum to 1, got {sum}");
}

wasm/crates/pf-wasm/src/pf.rs, trimmed. All seven tests pass.

Performance, measured rather than assumed

pnpm build:wasm pf-wasm produces a 40.9 kB binary, comfortably inside the workspace’s ~200 kB-per-crate guideline, and about as small as a crate with a PRNG, a 4-D EKF and a full particle cloud gets. The plan for this post predicted a 3–10× speedup over plain JS on the resampling and log-pdf inner loop. I measured it, rather than assumed it: a from-scratch JS port of the identical algorithm (same RNG, same model, same log-sum-exp, same systematic resampling), timed against the WASM module on the same 40-point sweep, in Node rather than a browser to remove rendering from the comparison:

Sweep	WASM	Plain JS	Speedup
$N_s$ sweep, 100→10 000, runs = 20	29.9 s	49.4 s	1.65×
Offset sweep, $N_s=2000$ , runs = 50	28.8 s	47.7 s	1.66×

Not 3–10×

The measured speedup is about 1.65×, not the 3–10× the plan predicted from post 19’s rendering numbers. I also tried the workspace’s opt-level = 3 override (the one imaging-wasm and segment-wasm use, temporarily, to measure it, reverted before publishing, since I don’t own wasm/Cargo.toml) and it made no measurable difference here. My best explanation: V8’s JIT is simply very good at this specific shape of code (tight scalar and typed-array loops with no allocation, no object headers, no bounds-check-defeating indirection), which is exactly the style this site’s JS convention already enforces (post 19’s particle cloud is one flat Float64Array, not an array of objects, for the same reason). WASM’s advantage tends to show up more on code that JS engines JIT poorly (branchy pixel kernels, SIMD-shaped work), not on this. Both numbers are still under a minute for the full budget, so the sweep is comfortably interactive either way; I built the WASM crate because the plan asked for it and because a ~1.6× win at this scale is real, but the honest lesson is that the “3–10×” framing came from a different post’s numbers, and this one’s own measurement doesn’t repeat it.

I did not need the cheap fallback (a precomputed sweep shipped as JSON, scrubbed rather than run): the live version comes in under a minute at the plan’s own honest budget, on both implementations, so the browser really does Monte Carlo it from a button press.

What I’d say now

§8’s conclusion, and I’d sign every sentence of it: the Kalman filter, in its extended form, is a genuinely useful numerical method that assumes the posterior is Gaussian, an assumption this whole series has spent four posts demolishing under the right conditions and vindicating under others. The particle filter approximates the posterior with samples instead, which works extremely well precisely where the Gaussian assumption fails, at the cost of more computation and no guarantee that any particle ever lands somewhere useful.

What I’d add, after building the sweep rather than reading about it: “how many particles” and “how wrong can the guess be” are not two separate questions with two separate answers. They’re the same question: how much does the true state actually cost to find, given this sensor, asked along two different axes. A cheap, well-observed problem (dataset A, or dataset C once the range measurement kills the ambiguity) needs neither many particles nor a good guess: the Kalman filter’s single Gaussian finds it in one matrix inversion. An expensive, badly-observed one (dataset B, especially along a velocity axis the sensor barely sees) needs both, and even then, as this post’s own $\dot x$ sweep showed, more particles buys you a lower probability of the cloud locking onto the wrong story, never a guarantee. The report’s title for its final section could as well have been the title of this whole series: match the filter, and its budget, to the actual shape of what you don’t know.