diff --git a/.github/workflows/test-coverage.yml b/.github/workflows/test-coverage.yml
new file mode 100644
index 0000000..dc78e9f
--- /dev/null
+++ b/.github/workflows/test-coverage.yml
@@ -0,0 +1,146 @@
+name: Test Coverage and Quality
+
+on:
+  push:
+    branches: [master, testing/infrastructure, testing/kalman-core, testing/kalman-advanced, testing/stats, testing/common, testing/other-modules]
+  pull_request:
+    branches: [master, testing/infrastructure, testing/kalman-core, testing/kalman-advanced, testing/stats, testing/common, testing/other-modules]
+
+jobs:
+  test:
+    name: Test on Python ${{ matrix.python-version }}
+    runs-on: ubuntu-latest
+    strategy:
+      matrix:
+        python-version: ['3.11', '3.12']
+
+    steps:
+      - uses: actions/checkout@v4
+
+      - name: Set up Python ${{ matrix.python-version }}
+        uses: actions/setup-python@v5
+        with:
+          python-version: ${{ matrix.python-version }}
+
+      - name: Install dependencies with uv
+        run: |
+          pip install uv
+          uv pip install -e ".[dev]" --system
+
+      - name: Run tests with coverage
+        run: |
+          pytest \
+            --cov=bayesian_filters \
+            --cov-report=xml \
+            --cov-report=html \
+            --cov-report=term-missing \
+            -v \
+            --tb=short
+
+      - name: Upload coverage reports to Codecov
+        uses: codecov/codecov-action@v4
+        with:
+          file: ./coverage.xml
+          flags: unittests
+          name: codecov-umbrella-py${{ matrix.python-version }}
+          fail_ci_if_error: false
+
+      - name: Archive coverage reports
+        if: always()
+        uses: actions/upload-artifact@v4
+        with:
+          name: coverage-report-py${{ matrix.python-version }}
+          path: htmlcov/
+
+  benchmarks:
+    name: Performance Benchmarks
+    runs-on: ubuntu-latest
+    if: github.event_name == 'push' && github.ref == 'refs/heads/master'
+
+    steps:
+      - uses: actions/checkout@v4
+
+      - name: Set up Python
+        uses: actions/setup-python@v5
+        with:
+          python-version: '3.11'
+
+      - name: Install dependencies with uv
+        run: |
+          pip install uv
+          uv pip install -e ".[dev]" --system
+
+      - name: Run benchmarks
+        run: |
+          pytest \
+            -m benchmark \
+            --benchmark-only \
+            --benchmark-json=.benchmarks/output.json \
+            -v
+
+      - name: Store benchmark result
+        uses: benchmark-action/github-action-benchmark@v1
+        with:
+          tool: 'pytest'
+          output-file-path: .benchmarks/output.json
+          github-token: ${{ secrets.GITHUB_TOKEN }}
+          auto-push: true
+        continue-on-error: true
+
+  quality:
+    name: Code Quality Checks
+    runs-on: ubuntu-latest
+
+    steps:
+      - uses: actions/checkout@v4
+
+      - name: Set up Python
+        uses: actions/setup-python@v5
+        with:
+          python-version: '3.11'
+
+      - name: Install dependencies with uv
+        run: |
+          pip install uv
+          uv pip install -e ".[dev]" --system
+
+      - name: Run type checking with pyright
+        run: pyright
+        continue-on-error: true
+
+      - name: Check code style with ruff
+        run: ruff check bayesian_filters tests --select E,W,F
+        continue-on-error: true
+
+  coverage-report:
+    name: Coverage Summary
+    runs-on: ubuntu-latest
+    needs: test
+    if: always()
+
+    steps:
+      - uses: actions/checkout@v4
+
+      - name: Set up Python
+        uses: actions/setup-python@v5
+        with:
+          python-version: '3.11'
+
+      - name: Install dependencies with uv
+        run: |
+          pip install uv
+          uv pip install -e ".[dev]" --system
+
+      - name: Generate coverage report
+        run: |
+          pytest \
+            --cov=bayesian_filters \
+            --cov-report=term-missing \
+            --cov-report=json \
+            -q
+
+      - name: Comment coverage on PR
+        if: github.event_name == 'pull_request'
+        uses: py-cov-action/python-coverage-comment-action@v3
+        with:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
diff --git a/TESTING.md b/TESTING.md
new file mode 100644
index 0000000..ef12f2c
--- /dev/null
+++ b/TESTING.md
@@ -0,0 +1,404 @@
+# Testing Guide for Bayesian Filters
+
+This document provides guidance on running tests, using testing utilities, and interpreting test results.
+
+## Quick Start
+
+### Running All Tests
+
+```bash
+pytest
+```
+
+### Running Tests with Coverage
+
+```bash
+pytest --cov=bayesian_filters --cov-report=html
+```
+
+This generates an HTML coverage report in `htmlcov/index.html`.
+
+### Running Specific Test Categories
+
+```bash
+# Unit tests only
+pytest -m unit
+
+# Integration tests
+pytest -m integration
+
+# Property-based tests (Hypothesis)
+pytest -m property
+
+# Numerical accuracy tests
+pytest -m numerical
+
+# Benchmarks
+pytest -m benchmark --benchmark-only
+
+# Slow tests (>1 second)
+pytest -m slow
+
+# Everything except slow tests
+pytest -m "not slow"
+```
+
+## Test Organization
+
+Tests are organized by module and test type:
+
+```
+bayesian_filters/
+├── kalman/
+│   ├── test_*.py              # Unit tests
+│   ├── test_*_property.py     # Property-based tests
+│   └── test_*_benchmark.py    # Performance benchmarks
+├── extended_kalman/
+├── unscented_kalman/
+├── common/
+└── stats/
+```
+
+## Using Testing Utilities
+
+### Fixtures
+
+Pre-configured filter fixtures are available via pytest:
+
+```python
+def test_with_fixture(kf_1d_fixture):
+    """Test using a 1D Kalman filter fixture."""
+    kf = kf_1d_fixture
+    # kf has state [position, velocity], dt=0.1
+    assert kf.x.shape == (2, 1)
+
+def test_2d_filter(kf_2d_fixture):
+    """Test using a 2D Kalman filter fixture."""
+    kf = kf_2d_fixture
+    # kf has state [x, vx, y, vy]
+    assert kf.x.shape == (4, 1)
+```
+
+**Available fixtures:**
+- `kf_1d_fixture`: 1D constant-velocity Kalman filter
+- `kf_2d_fixture`: 2D constant-velocity Kalman filter
+- `ekf_2d_fixture`: 2D Extended Kalman filter with range measurement
+- `ukf_3d_fixture`: 3D Unscented Kalman filter
+- `sensor_sim_fixture`: Simple 1D sensor simulator
+
+### Helper Functions
+
+Test assertion utilities for validating filter properties:
+
+```python
+from bayesian_filters.testing_utils import (
+    assert_matrix_psd,
+    assert_matrix_symmetric,
+    assert_filter_stable,
+)
+
+def test_covariance_properties(kf_1d_fixture):
+    """Verify filter covariance matrix properties."""
+    kf = kf_1d_fixture
+
+    # Check covariance is positive semi-definite
+    assert_matrix_psd(kf.P, name="Initial covariance")
+
+    # Check covariance is symmetric
+    assert_matrix_symmetric(kf.P, name="Initial covariance")
+
+    # Check filter stability
+    kf.predict()
+    assert_filter_stable(kf.P, max_variance=1e10)
+```
+
+### Numerical Solutions
+
+Analytical solutions for validation:
+
+```python
+from bayesian_filters.testing_utils import ConstantVelocitySolution
+
+def test_against_analytical_solution():
+    """Compare filter estimate against analytical solution."""
+    dt = 0.1
+    x0, v0 = 0.0, 1.0
+    t = 1.0
+
+    # Analytical position at t=1.0
+    analytical_pos = ConstantVelocitySolution.position_at_time(t, x0, v0)
+
+    # Compare with filter estimate
+    assert abs(filter_estimate - analytical_pos) < 0.1
+```
+
+## Test Types
+
+### Unit Tests
+
+Fast, isolated tests for individual functions or methods.
+
+```python
+@pytest.mark.unit
+def test_kalman_gain_computation(kf_1d_fixture):
+    """Test Kalman gain computation."""
+    kf = kf_1d_fixture
+    # Test implementation
+    assert K.shape == (2, 1)
+```
+
+### Property-Based Tests
+
+Use Hypothesis to generate random test cases and verify invariants:
+
+```python
+from hypothesis import given, strategies as st
+
+@pytest.mark.property
+@given(
+    measurement=st.floats(min_value=-100, max_value=100),
+    noise_cov=st.floats(min_value=0.1, max_value=10.0),
+)
+def test_filter_convergence(measurement, noise_cov):
+    """Test that filter estimates converge with Hypothesis."""
+    kf = KalmanFilter(dim_x=2, dim_z=1)
+    # ... setup ...
+
+    # Filter should produce reasonable estimates
+    kf.update(measurement)
+    assert np.isfinite(kf.x).all()
+    assert kf.P[0, 0] > 0
+```
+
+### Numerical Accuracy Tests
+
+Verify filter accuracy against known solutions:
+
+```python
+@pytest.mark.numerical
+def test_constant_velocity_tracking():
+    """Test filter tracks constant velocity motion accurately."""
+    # Known motion
+    ground_truth = ConstantVelocitySolution.position_at_time(t=1.0)
+
+    # Filter estimate
+    kf = KalmanFilter(dim_x=2, dim_z=1)
+    kf.x = np.array([[0.0], [1.0]])
+    kf.predict()
+
+    # Error should be small
+    error = abs(kf.x[0, 0] - ground_truth)
+    assert error < 0.01
+```
+
+### Performance Benchmarks
+
+Measure filter performance and detect regressions:
+
+```python
+@pytest.mark.benchmark
+def test_kalman_filter_predict_benchmark(benchmark, kf_1d_fixture):
+    """Benchmark Kalman filter prediction speed."""
+    kf = kf_1d_fixture
+
+    # Measure predict() speed
+    result = benchmark(kf.predict)
+
+    # Should complete quickly
+    assert benchmark.stats.mean < 0.001  # < 1 ms
+```
+
+## Coverage Requirements
+
+Test coverage targets by module:
+
+| Module | Target | Notes |
+|--------|--------|-------|
+| kalman.KalmanFilter | 95% | Core algorithm |
+| kalman.ExtendedKalmanFilter | 90% | Includes nonlinearity |
+| kalman.UnscentedKalmanFilter | 90% | Sigma point algorithm |
+| common | 85% | Utility functions |
+| stats | 80% | Statistical functions |
+| Other modules | 75% | Supporting code |
+
+## Running Coverage Analysis
+
+### Generate Coverage Report
+
+```bash
+pytest --cov=bayesian_filters --cov-report=html
+open htmlcov/index.html
+```
+
+### Show Missing Lines
+
+```bash
+pytest --cov=bayesian_filters --cov-report=term-missing
+```
+
+### Coverage for Specific Module
+
+```bash
+pytest --cov=bayesian_filters.kalman --cov-report=term-missing
+```
+
+## Parallel Test Execution
+
+Run tests in parallel for faster feedback:
+
+```bash
+# Run on all CPU cores
+pytest -n auto
+
+# Run on specific number of cores
+pytest -n 4
+```
+
+## Continuous Integration
+
+Tests are automatically run on GitHub Actions for:
+- All push events to testing branches
+- All pull requests
+- Multiple Python versions (3.9, 3.10, 3.11, 3.12)
+- Coverage collection and reporting
+- Performance benchmarks on master branch
+
+See `.github/workflows/test-coverage.yml` for full CI/CD configuration.
+
+## Debugging Test Failures
+
+### Verbose Output
+
+```bash
+pytest -vv --tb=long test_file.py::test_function
+```
+
+### Drop into Debugger
+
+```bash
+pytest --pdb test_file.py::test_function
+```
+
+### Show Print Statements
+
+```bash
+pytest -s test_file.py::test_function
+```
+
+### Run with Random Seed
+
+```bash
+# Reproducible randomization
+pytest --seed=12345
+
+# Get seed from last failure
+pytest --lastfailed --seed=12345
+```
+
+## Writing New Tests
+
+### Test Structure
+
+```python
+"""Tests for kalman.py module.
+
+This module tests the KalmanFilter class including:
+- Filter initialization
+- Prediction step
+- Update step
+- Covariance matrix properties
+"""
+
+import pytest
+import numpy as np
+from bayesian_filters.kalman import KalmanFilter
+from bayesian_filters.testing_utils import assert_matrix_psd
+
+@pytest.mark.unit
+class TestKalmanFilterInitialization:
+    """Test KalmanFilter initialization."""
+
+    def test_state_initialization(self, kf_1d_fixture):
+        """Verify initial state is set correctly."""
+        kf = kf_1d_fixture
+        assert kf.x.shape == (2, 1)
+        assert kf.x[0, 0] == 0.0
+        assert kf.x[1, 0] == 1.0
+
+    def test_covariance_initialization(self, kf_1d_fixture):
+        """Verify initial covariance is positive definite."""
+        kf = kf_1d_fixture
+        assert_matrix_psd(kf.P)
+
+@pytest.mark.unit
+def test_kalman_filter_predict(kf_1d_fixture):
+    """Test prediction step maintains covariance properties."""
+    kf = kf_1d_fixture
+
+    # Predict
+    kf.predict()
+
+    # Covariance should still be valid
+    assert_matrix_psd(kf.P)
+    assert np.isfinite(kf.x).all()
+    assert np.isfinite(kf.P).all()
+```
+
+### Guidelines
+
+1. **Use fixtures** for common setup (kf_1d_fixture, sensor_sim_fixture, etc.)
+2. **Test invariants** - verify properties that should always hold
+3. **Use markers** - @pytest.mark.unit, @pytest.mark.numerical, etc.
+4. **Clear names** - test function names should describe what they test
+5. **Assertions** - use helper functions from testing_utils when possible
+6. **Avoid randomness** - use fixed seeds or Hypothesis for controlled randomness
+
+## Interpreting Coverage Reports
+
+The HTML coverage report (`htmlcov/index.html`) shows:
+
+- **Green lines**: Covered by tests
+- **Red lines**: Not covered by tests
+- **Yellow lines**: Partially covered
+- **Missing branches**: Shows which if/else branches aren't tested
+
+Aim for:
+- **High line coverage** (80-95%)
+- **Branch coverage** for critical paths (>90%)
+- **All public APIs** tested
+
+## Common Issues
+
+### "No tests ran"
+```bash
+# Check test discovery
+pytest --collect-only
+
+# Verify test file naming (test_*.py or *_test.py)
+# Verify test function naming (test_* prefix)
+```
+
+### "Fixtures not found"
+```bash
+# Check conftest.py is in project root
+# Verify fixture is exported in __init__.py
+# Use pytest --fixtures to list available fixtures
+pytest --fixtures | grep kalman
+```
+
+### "Coverage missing"
+```bash
+# Reinstall in development mode
+pip install -e .[dev]
+
+# Check pyproject.toml has omit patterns correct
+# Run with --no-cov-on-fail
+pytest --cov --no-cov-on-fail
+```
+
+## Resources
+
+- [pytest documentation](https://docs.pytest.org/)
+- [Hypothesis documentation](https://hypothesis.readthedocs.io/)
+- [pytest-benchmark](https://pytest-benchmark.readthedocs.io/)
+- Testing Plan: See `TESTING_PLAN.md` for comprehensive testing strategy
diff --git a/TESTING_PLAN.md b/TESTING_PLAN.md
new file mode 100644
index 0000000..d0caa91
--- /dev/null
+++ b/TESTING_PLAN.md
@@ -0,0 +1,480 @@
+# Comprehensive Unit Testing Plan - Worktree/Branch Structure
+
+**Status:** Approved and In Progress
+**Created:** 2025-10-25
+**Target Coverage:** 80-90% line coverage
+**Total Estimated New Tests:** 400-550 tests
+
+## Overview
+
+Create a systematic testing infrastructure organized by module, targeting 80-90% line coverage with emphasis on property-based testing, numerical accuracy, parametrization, and performance benchmarks.
+
+---
+
+## Branch/Worktree Structure
+
+### 1. testing/infrastructure (Foundation - Start Here)
+
+**Worktree:** `/tmp/bayesian-filters-testing-infrastructure`
+**Base Branch:** `master`
+**Status:** IN PROGRESS
+
+#### Scope
+- Add pytest-cov for coverage reporting
+- Add hypothesis for property-based testing
+- Add pytest-benchmark for performance tests
+- Configure pyproject.toml with pytest settings
+- Add pytest-xdist for parallel test execution
+- Create GitHub Actions workflow for test coverage reporting
+- Add test utilities and fixtures in `bayesian_filters/testing_utils/`
+- Document testing standards and conventions
+
+#### Deliverables
+- `pyproject.toml` updated with `[tool.pytest.ini_options]` and `[tool.coverage.run]`
+- `.github/workflows/test-coverage.yml` for CI coverage reporting
+- `bayesian_filters/testing_utils/__init__.py`
+- `bayesian_filters/testing_utils/fixtures.py` with common test fixtures
+- `bayesian_filters/testing_utils/numerical.py` with analytical solutions for comparison
+- `bayesian_filters/testing_utils/helpers.py` with test utilities
+- `TESTING.md` documentation in root
+- `conftest.py` at root level for pytest plugins
+
+#### Key Dependencies
+```
+pytest>=7.0
+pytest-cov>=4.0
+hypothesis>=6.0
+pytest-benchmark>=4.0
+pytest-xdist>=3.0
+pytest-mpl>=0.15 (for visual regression)
+```
+
+#### Implementation Checklist
+- [ ] Update pyproject.toml with pytest configuration
+- [ ] Create testing_utils package structure
+- [ ] Implement basic fixtures (kf_1d, ekf_2d, ukf_3d, etc.)
+- [ ] Implement analytical solutions library
+- [ ] Create conftest.py with pytest plugins
+- [ ] Add GitHub Actions test-coverage.yml workflow
+- [ ] Write TESTING.md documentation
+- [ ] Test infrastructure setup
+
+---
+
+### 2. testing/kalman-core (PRIORITY - Core Filters)
+
+**Worktree:** `/tmp/bayesian-filters-testing-kalman-core`
+**Base Branch:** `testing/infrastructure`
+**Status:** PENDING
+
+#### Modules to Test
+- **kalman_filter.py** (KalmanFilter class) - CRITICAL
+- **EKF.py** (ExtendedKalmanFilter) - CRITICAL
+- **UKF.py** (UnscentedKalmanFilter) - CRITICAL
+
+#### Current Test State
+| Module | File | Tests | Lines | Status |
+|--------|------|-------|-------|--------|
+| KalmanFilter | test_kf.py | 14 | 738 | Decent |
+| EKF | test_ekf.py | 1 | 122 | MINIMAL |
+| UKF | test_ukf.py | 16 | 1108 | Good |
+
+#### Testing Focus
+
+##### 1. Unit Tests
+- Test each method independently (predict, update, batch_filter, etc.)
+- Edge cases: singular matrices, zero covariance, dimension mismatches
+- Error conditions: invalid inputs, numerical instability
+- State shape handling (row vs column vectors)
+- Matrix dimension validation
+
+##### 2. Property-Based Tests (Hypothesis)
+- Generate random valid filter configurations
+- Verify conservation properties (covariance remains positive semi-definite)
+- Test invariants (symmetry of P, positive definiteness after update)
+- Commutative properties where applicable
+
+##### 3. Numerical Accuracy Tests
+- Compare against analytical solutions (constant velocity, constant acceleration)
+- Verify against published benchmark results
+- Test numerical stability with ill-conditioned matrices
+- Cross-validation between different filter implementations
+
+##### 4. Parametrized Tests
+- Different state dimensions (1D, 2D, 3D, 4D, 10D)
+- Different measurement dimensions
+- Various noise levels (high, medium, low SNR)
+- Different control signal scenarios
+- Batch vs sequential filtering equivalence
+
+##### 5. Performance Benchmarks
+- Baseline filter update/predict times
+- Scaling with state dimension (O(n²) vs O(n³))
+- Memory usage profiles
+- Vectorization effectiveness
+
+#### Estimated New Tests
+- KalmanFilter: 50-60 new tests
+- EKF: 40-50 new tests
+- UKF: 30-40 new tests
+- **Total: 120-150 new tests**
+
+#### New Test File Targets
+- `bayesian_filters/kalman/tests/test_kf_comprehensive.py` - Extended KF tests
+- `bayesian_filters/kalman/tests/test_ekf_comprehensive.py` - Extended EKF tests
+- `bayesian_filters/kalman/tests/test_ukf_comprehensive.py` - Extended UKF tests
+
+---
+
+### 3. testing/kalman-advanced (Advanced Filters)
+
+**Worktree:** `/tmp/bayesian-filters-testing-kalman-advanced`
+**Base Branch:** `testing/infrastructure`
+**Status:** PENDING
+
+#### Modules to Test
+- **CubatureKalmanFilter.py** (test_ckf.py: 1 test → expand to 20-30)
+- **ensemble_kalman_filter.py** (test_enkf.py: 2 tests → expand to 15-25)
+- **IMM.py** (test_imm.py: 2 tests → expand to 20-30)
+- **information_filter.py** (test_information.py: 12 tests → expand to 25-35)
+- **square_root.py** (test_sqrtkf.py: 1 test → expand to 20-30)
+- **fading_memory.py** (test_fm.py: 1 test → expand to 15-25)
+- **fixed_lag_smoother.py** (test_fls.py: 2 tests → expand to 15-25)
+- **mmae.py** (test_mmae.py: 1 test → expand to 20-30)
+- **RTS Smoother** (test_rts.py: 1 test → expand to 20-30)
+- **Sensor Fusion** (test_sensor_fusion.py: 1 test → expand to 15-25)
+
+#### Testing Approach
+- Same comprehensive approach as core filters
+- Cross-validation between filter types (same problem, different filters)
+- Comparison with standard KF/EKF/UKF where applicable
+- Filter-specific edge cases:
+  - EnKF: ensemble size, covariance localization
+  - IMM: model switching, mixing weights
+  - CKF: cubature point distribution
+  - Information Filter: information matrix singularity
+  - RTS: backward pass correctness
+  - Fading Memory: forgetting factor effects
+
+#### Estimated New Tests
+- Per module: 20-30 new tests
+- **Total: 160-240 new tests**
+
+---
+
+### 4. testing/stats (Statistical Functions)
+
+**Worktree:** `/tmp/bayesian-filters-testing-stats`
+**Base Branch:** `testing/infrastructure`
+**Status:** PENDING
+
+#### Module
+- **stats/stats.py**
+- Current: test_stats.py: 5 tests, 301 lines
+
+#### Functions Needing Thorough Testing
+- Gaussian operations (gaussian, multivariate_gaussian, mul, add, mul_pdf)
+- Mahalanobis distance
+- Log-likelihood and likelihood
+- Plotting functions (plot_gaussian_pdf, plot_covariance, plot_gaussian_cdf, etc.)
+- Statistical utilities (NEES, norm_cdf)
+- Covariance ellipse calculations
+
+#### Testing Focus
+- Numerical accuracy against scipy.stats
+- Edge cases:
+  - Singular covariances
+  - High dimensions
+  - Extreme values (very small/large covariances)
+  - Near-zero variances
+- Property-based tests for Gaussian algebra:
+  - Multiplication is associative where defined
+  - Addition is commutative
+  - Mahalanobis distance properties
+- Visual regression testing for plotting functions (using pytest-mpl)
+- Comparison with scipy implementations
+
+#### Estimated New Tests
+- **30-40 new tests**
+
+#### Coverage Targets
+- 95%+ coverage (mostly plotting, validation code)
+
+---
+
+### 5. testing/common (Utility Functions)
+
+**Worktree:** `/tmp/bayesian-filters-testing-common`
+**Base Branch:** `testing/infrastructure`
+**Status:** PENDING
+
+#### Modules
+- **helpers.py** (test_helpers.py: 7 tests → expand to 20-30)
+- **discretization.py** (test_discretization.py: 3 tests → expand to 15-20)
+- **kinematic.py** (test_kinematic.py: 29 test methods → expand to 40-50)
+
+#### Testing Focus
+- Matrix operations accuracy
+- Q_discrete_white_noise correctness (compare to analytical Wiener process)
+- Kinematic model validation against analytical solutions
+- Helper function edge cases
+- Block diagonal and other matrix utilities
+
+#### Specific Tests
+- Discretization accuracy for different dt values
+- Kinematic state transition equivalence
+- Saver class functionality
+- Q matrix positive definiteness
+
+#### Estimated New Tests
+- **20-30 new tests**
+
+---
+
+### 6. testing/other-modules (Remaining Modules)
+
+**Worktree:** `/tmp/bayesian-filters-testing-other`
+**Base Branch:** `testing/infrastructure`
+**Status:** PENDING
+
+#### Modules
+- **discrete_bayes/** (test_discrete_bayes.py: 1 test → expand to 10-15)
+- **gh/** (test_gh.py: 4 tests → expand to 15-20)
+- **hinfinity/** (test_hinfinity.py: 1 test → expand to 10-15)
+- **leastsq/** (test_lsq.py: 6 tests → expand to 20-25)
+- **memory/** (test_fading_memory.py: 1 test → expand to 10-15)
+- **monte_carlo/** (test_resampling.py: 23 test methods → expand to 50-60)
+
+#### Testing Focus
+- Module-specific algorithms
+- Cross-module integration
+- Edge cases and error handling
+- Algorithm correctness against publications
+
+#### Estimated New Tests
+- **30-50 new tests**
+
+---
+
+## Testing Standards
+
+### Test Organization
+
+Each test file should contain sections in this order:
+```python
+# 1. Imports and setup
+# 2. Fixtures and helpers
+# 3. Unit tests (test_<method_name>_<scenario>)
+# 4. Integration tests (test_integration_<workflow>)
+# 5. Property-based tests (test_property_<invariant>)
+# 6. Parametrized tests (using @pytest.mark.parametrize)
+# 7. Performance benchmarks (benchmark_<operation>)
+```
+
+### Coverage Requirements
+- **Target:** 80-90% line coverage per module
+- **Critical paths:** 100% coverage for predict/update methods
+- **Error handling:** All error conditions must be tested
+- **Edge cases:** Documented and tested
+- **Numerical stability:** Tested with condition numbers >1e6
+
+### Naming Conventions
+```python
+test_<method_name>_<scenario>              # Unit tests
+    e.g., test_predict_constant_velocity
+    e.g., test_update_dimension_mismatch
+
+test_integration_<workflow_name>           # Integration tests
+    e.g., test_integration_batch_filtering
+    e.g., test_integration_multi_model_switching
+
+test_property_<invariant_name>             # Property-based (Hypothesis)
+    e.g., test_property_covariance_remains_psd
+    e.g., test_property_symmetry_preserved
+
+benchmark_<operation>                      # Performance tests
+    e.g., benchmark_predict_1d_filter
+    e.g., benchmark_update_scaling_with_dimension
+```
+
+### Test Markers
+```python
+@pytest.mark.unit           # Fast, isolated tests
+@pytest.mark.integration    # Tests multiple components
+@pytest.mark.slow          # Tests that take >1 second
+@pytest.mark.benchmark     # Performance benchmarks
+@pytest.mark.parametrize   # Multiple parameter sets
+```
+
+### Assertion Style
+- Use `pytest.approx()` for float comparisons (not `==`)
+- Use `np.allclose()` for array comparisons
+- Document tolerance in comments
+- Use descriptive assertion messages
+
+---
+
+## Implementation Timeline
+
+### Week 1: Infrastructure
+- [ ] Create testing/infrastructure branch
+- [ ] Update pyproject.toml
+- [ ] Create testing_utils package
+- [ ] Implement fixtures and helpers
+- [ ] Set up GitHub Actions workflow
+- [ ] Write TESTING.md documentation
+
+### Week 2-3: Core Kalman Filters
+- [ ] Expand test_kf.py to 50-60 tests
+- [ ] Expand test_ekf.py from 1 to 40-50 tests
+- [ ] Expand test_ukf.py to 30-40 tests
+- [ ] Add property-based tests
+- [ ] Add benchmarks
+- [ ] Achieve 90%+ coverage
+
+### Week 4: Advanced Kalman Filters
+- [ ] Expand all advanced filter tests
+- [ ] Cross-validation tests
+- [ ] Algorithm correctness verification
+- [ ] Achieve 85%+ coverage
+
+### Week 5: Stats & Common
+- [ ] Expand stats/stats.py tests
+- [ ] Expand common utility tests
+- [ ] Add numerical accuracy tests
+- [ ] Achieve 90%+ coverage
+
+### Week 6: Other Modules & Integration
+- [ ] Expand remaining module tests
+- [ ] Add integration tests
+- [ ] Full system validation
+- [ ] Coverage reporting
+
+---
+
+## Merge Strategy
+
+Each testing branch will be reviewed and merged independently:
+
+1. **Create PR** from testing branch to master
+2. **Requirements:**
+   - 80%+ coverage in modified modules
+   - All tests pass
+   - Performance benchmarks must not regress >10%
+   - Code review approval
+3. **Merge** when all requirements met
+4. **Parallel Development:** Branches can be worked on in parallel
+
+### PR Template
+```markdown
+## Testing Enhancement for [Module]
+
+### Coverage
+- Current: X%
+- After: Y%
+- New tests: N
+
+### Test Types
+- [ ] Unit tests
+- [ ] Property-based tests
+- [ ] Integration tests
+- [ ] Numerical accuracy tests
+- [ ] Performance benchmarks
+
+### Validation
+- [ ] All tests pass
+- [ ] Coverage increased
+- [ ] No performance regressions
+- [ ] Documentation updated
+```
+
+---
+
+## Performance Baseline
+
+Benchmarks will establish baselines for:
+- KalmanFilter.predict() - 1D to 10D states
+- KalmanFilter.update() - 1D to 10D measurements
+- EKF operations with nonlinear functions
+- UKF operations with different sigma points
+- Batch filtering performance
+
+Regression threshold: 10% slowdown triggers investigation
+
+---
+
+## Testing Utilities Overview
+
+### `testing_utils/fixtures.py`
+Common pytest fixtures:
+- `kf_1d` - 1D Kalman filter
+- `kf_2d` - 2D Kalman filter
+- `ekf_2d` - 2D Extended Kalman filter
+- `ukf_3d` - 3D Unscented Kalman filter
+- `sensor_sim` - Simple sensor simulator
+- `noisy_measurements` - Generate measurement sequences
+
+### `testing_utils/numerical.py`
+Analytical solutions and comparisons:
+- Constant velocity filter analytical solution
+- Constant acceleration filter solution
+- Linear system steady state
+- Benchmark filter configurations
+- Published test cases
+
+### `testing_utils/helpers.py`
+Test utilities:
+- Matrix verification functions
+- Filter validation helpers
+- Performance measurement utilities
+- Comparison against scipy/numpy
+
+---
+
+## Success Criteria
+
+✅ **Complete when:**
+1. All 6 branches created and merged
+2. Total coverage: 80-90% across all modules
+3. Critical paths (predict/update): 100%
+4. 400-550 new tests implemented
+5. All property-based tests passing
+6. Benchmarks established and no regressions
+7. TESTING.md updated with results
+8. CI/CD pipeline reporting coverage
+9. Documentation complete
+
+---
+
+## References & Resources
+
+### Testing Frameworks
+- [pytest Documentation](https://docs.pytest.org/)
+- [Hypothesis Guide](https://hypothesis.readthedocs.io/)
+- [pytest-benchmark](https://pytest-benchmark.readthedocs.io/)
+
+### Kalman Filter Testing
+- Welch & Bishop: "An Introduction to the Kalman Filter"
+- Published benchmark datasets
+- Cross-reference with reference implementations
+
+### Numerical Testing
+- IEEE 754 floating point standards
+- Condition number analysis
+- Matrix stability considerations
+
+---
+
+## Notes
+
+- Tests should be deterministic (use fixed seeds for randomness)
+- Avoid hardcoded magic numbers - use descriptive constants
+- Document assumptions and edge cases in docstrings
+- Use type hints in test functions where helpful
+- Keep tests focused and independent
+
+---
+
+**Last Updated:** 2025-10-25
+**Plan Version:** 1.0
+**Approved:** Yes
diff --git a/bayesian_filters/kalman/tests/test_kf_comprehensive.py b/bayesian_filters/kalman/tests/test_kf_comprehensive.py
new file mode 100644
index 0000000..ba9bb63
--- /dev/null
+++ b/bayesian_filters/kalman/tests/test_kf_comprehensive.py
@@ -0,0 +1,701 @@
+"""Comprehensive tests for KalmanFilter class.
+
+This module provides extensive testing of the KalmanFilter implementation
+across four categories:
+
+1. Analytical Solution Tests: Verify against known analytical solutions
+2. Property-Based Invariant Tests: Use Hypothesis to verify invariants
+3. Edge Case and Error Handling: Test boundary conditions
+4. Numerical Stability: Test long-running and ill-conditioned scenarios
+
+These tests complement the existing test_kf.py by adding coverage for:
+- Numerical accuracy against known solutions
+- Property-based testing of matrix invariants
+- Comprehensive edge case handling
+- Long-term stability and convergence
+"""
+
+import numpy as np
+import pytest
+from hypothesis import given, strategies as st, settings, HealthCheck
+from bayesian_filters.kalman import KalmanFilter
+from bayesian_filters.common import Q_discrete_white_noise
+from bayesian_filters.testing_utils import (
+    assert_matrix_psd,
+    assert_matrix_symmetric,
+    assert_filter_stable,
+    ConstantVelocitySolution,
+    ConstantAccelerationSolution,
+)
+
+
+# ============================================================================
+# PART 1: ANALYTICAL SOLUTION TESTS
+# ============================================================================
+
+
+class TestKalmanFilterAnalyticalSolutions:
+    """Test KalmanFilter against known analytical solutions."""
+
+    @pytest.mark.numerical
+    def test_constant_velocity_perfect_measurement(self):
+        """Test tracking with perfect measurement (R=0, small Q).
+
+        Setup: 1D constant velocity system with R=0 (perfect measurement), Q small.
+        Expected: Position estimate should match measurements exactly.
+        """
+        dt = 0.1
+        kf = KalmanFilter(dim_x=2, dim_z=1)
+
+        # Setup constant velocity model
+        kf.x = np.array([[0.0], [1.0]])
+        kf.F = np.array([[1.0, dt], [0.0, 1.0]])
+        kf.H = np.array([[1.0, 0.0]])
+        kf.P = np.eye(2) * 100.0
+        kf.Q = np.eye(2) * 1e-8  # Very small process noise
+        kf.R = np.array([[1e-10]])  # Very small measurement noise
+
+        # Simulate constant velocity motion
+        for t in range(10):
+            # Analytical position
+            z_true = ConstantVelocitySolution.position_at_time(t * dt, x0=0.0, v0=1.0)
+
+            # Predict and update with near-perfect measurement
+            kf.predict()
+            kf.update(z_true)
+
+            # Position should match measurements closely
+            assert np.isclose(kf.x[0, 0], z_true, atol=0.01)
+
+    @pytest.mark.numerical
+    def test_constant_velocity_with_noise(self):
+        """Test convergence with process and measurement noise.
+
+        Setup: 1D constant velocity system with Q>0, R>0.
+        Expected: Filter state should track within 2-3 sigma of truth.
+        """
+        dt = 0.1
+        kf = KalmanFilter(dim_x=2, dim_z=1)
+
+        # Setup constant velocity model with noise
+        kf.x = np.array([[0.0], [1.0]])
+        kf.F = np.array([[1.0, dt], [0.0, 1.0]])
+        kf.H = np.array([[1.0, 0.0]])
+        kf.P = np.eye(2) * 100.0
+        kf.Q = Q_discrete_white_noise(dim=2, dt=dt, var=0.01)
+        kf.R = np.array([[1.0]])
+
+        # Seed for reproducibility
+        np.random.seed(42)
+
+        # Simulate with noise
+        errors = []
+        for t in range(50):
+            # Ground truth
+            z_true = ConstantVelocitySolution.position_at_time(t * dt, x0=0.0, v0=1.0)
+
+            # Noisy measurement
+            z_measured = z_true + np.random.normal(0, 1.0)
+
+            kf.predict()
+            kf.update(z_measured)
+
+            error = np.abs(kf.x[0, 0] - z_true)
+            errors.append(error)
+
+        # Final error should converge
+        assert np.mean(errors[-10:]) < np.mean(errors[:10])
+
+    @pytest.mark.numerical
+    def test_constant_acceleration(self):
+        """Test constant acceleration model.
+
+        Setup: 1D constant acceleration system.
+        Expected: Position estimates should be reasonable.
+        """
+        dt = 0.1
+        kf = KalmanFilter(dim_x=3, dim_z=1)
+
+        # Setup constant acceleration model
+        kf.x = np.array([[0.0], [0.0], [1.0]])
+        kf.F = np.array([[1.0, dt, 0.5 * dt**2], [0.0, 1.0, dt], [0.0, 0.0, 1.0]])
+        kf.H = np.array([[1.0, 0.0, 0.0]])
+        kf.P = np.eye(3) * 100.0
+        kf.Q = Q_discrete_white_noise(dim=3, dt=dt, var=0.01)
+        kf.R = np.array([[0.1]])
+
+        np.random.seed(42)
+
+        for t in range(20):
+            # Analytical solution
+            z_true = ConstantAccelerationSolution.position_at_time(t * dt, x0=0.0, v0=0.0, a=1.0)
+
+            z_measured = z_true + np.random.normal(0, np.sqrt(0.1))
+
+            kf.predict()
+            kf.update(z_measured)
+
+            # Check estimates are reasonable (within measurement noise range)
+            assert np.abs(kf.x[0, 0] - z_true) < 1.0
+            assert np.isfinite(kf.x).all()
+
+    @pytest.mark.numerical
+    def test_2d_tracking(self):
+        """Test 2D constant velocity tracking.
+
+        Setup: 2D constant velocity system.
+        Expected: Both x and y coordinates should be tracked reasonably.
+        """
+        dt = 0.1
+        kf = KalmanFilter(dim_x=4, dim_z=2)
+
+        # Setup 2D constant velocity model
+        kf.x = np.array([[0.0], [1.0], [0.0], [1.0]])
+        F = np.eye(4)
+        F[0, 1] = dt
+        F[2, 3] = dt
+        kf.F = F
+
+        kf.H = np.array([[1.0, 0.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0]])
+        kf.P = np.eye(4) * 100.0
+        kf.Q = np.eye(4) * 0.01
+        kf.R = np.eye(2)
+
+        np.random.seed(42)
+
+        for t in range(10):
+            # Ground truth (moving diagonally)
+            z_x = t * dt
+            z_y = t * dt
+
+            z_measured = np.array([[z_x], [z_y]]) + np.random.randn(2, 1)
+
+            kf.predict()
+            kf.update(z_measured)
+
+            # Both coordinates should be tracked reasonably
+            assert np.abs(kf.x[0, 0] - z_x) < 1.5
+            assert np.abs(kf.x[2, 0] - z_y) < 1.5
+            assert np.isfinite(kf.x).all()
+
+    @pytest.mark.numerical
+    def test_steady_state_convergence(self):
+        """Test that filter covariance settles after initialization.
+
+        Setup: 1D constant velocity system with constant noise.
+        Expected: Covariance should settle (not growing indefinitely).
+        """
+        dt = 0.1
+        kf = KalmanFilter(dim_x=2, dim_z=1)
+
+        kf.x = np.array([[0.0], [1.0]])
+        kf.F = np.array([[1.0, dt], [0.0, 1.0]])
+        kf.H = np.array([[1.0, 0.0]])
+        kf.P = np.eye(2) * 100.0
+        kf.Q = Q_discrete_white_noise(dim=2, dt=dt, var=0.1)
+        kf.R = np.array([[1.0]])
+
+        # Run filter to convergence
+        covariances = []
+        for _ in range(100):
+            kf.predict()
+            kf.update(np.array([[0.0]]))
+            covariances.append(np.trace(kf.P))
+
+        # Later traces should be stable (bounded)
+        trace_early = np.mean(covariances[:10])
+        trace_late = np.mean(covariances[-10:])
+
+        # Both should be finite
+        assert np.isfinite(trace_early)
+        assert np.isfinite(trace_late)
+        # Later trace should not be growing rapidly
+        assert trace_late < trace_early * 10
+
+
+# ============================================================================
+# PART 2: PROPERTY-BASED INVARIANT TESTS
+# ============================================================================
+
+
+class TestKalmanFilterInvariants:
+    """Test invariants that must always hold for KalmanFilter."""
+
+    @pytest.mark.property
+    def test_covariance_remains_symmetric(self, kf_1d_fixture):
+        """Property: Covariance P must always remain symmetric."""
+        kf = kf_1d_fixture
+
+        for _ in range(20):
+            kf.predict()
+            assert_matrix_symmetric(kf.P, name="P after predict")
+
+            kf.update(np.array([[0.0]]))
+            assert_matrix_symmetric(kf.P, name="P after update")
+
+    @pytest.mark.property
+    def test_covariance_positive_semidefinite(self, kf_1d_fixture):
+        """Property: Covariance P must always remain positive semi-definite."""
+        kf = kf_1d_fixture
+
+        for _ in range(20):
+            kf.predict()
+            assert_matrix_psd(kf.P, name="P after predict")
+
+            kf.update(np.array([[np.random.randn()]]))
+            assert_matrix_psd(kf.P, name="P after update")
+
+    @pytest.mark.property
+    @given(
+        measurements=st.lists(
+            st.floats(min_value=-100, max_value=100, allow_nan=False),
+            min_size=1,
+            max_size=50,
+        )
+    )
+    @settings(suppress_health_check=[HealthCheck.function_scoped_fixture])
+    def test_covariance_reduces_after_update(self, measurements, kf_1d_fixture):
+        """Property: Trace of P should decrease after update (measurement reduces uncertainty)."""
+        kf = kf_1d_fixture
+
+        for z in measurements:
+            kf.predict()
+            trace_before_update = np.trace(kf.P)
+
+            kf.update(np.array([[z]]))
+            trace_after_update = np.trace(kf.P)
+
+            # Update should not increase uncertainty (within numerical tolerance)
+            assert trace_after_update <= trace_before_update + 1e-10
+
+    @pytest.mark.property
+    @given(
+        measurements=st.lists(
+            st.floats(min_value=-100, max_value=100, allow_nan=False),
+            min_size=1,
+            max_size=50,
+        )
+    )
+    @settings(suppress_health_check=[HealthCheck.function_scoped_fixture])
+    def test_state_remains_finite(self, measurements, kf_1d_fixture):
+        """Property: State vector x must always contain finite values."""
+        kf = kf_1d_fixture
+
+        for z in measurements:
+            kf.predict()
+            assert np.isfinite(kf.x).all(), "State contains NaN or Inf after predict"
+
+            kf.update(np.array([[z]]))
+            assert np.isfinite(kf.x).all(), "State contains NaN or Inf after update"
+
+    @pytest.mark.property
+    def test_covariance_remains_finite(self, kf_1d_fixture):
+        """Property: Covariance P must always contain finite values."""
+        kf = kf_1d_fixture
+
+        for _ in range(50):
+            kf.predict()
+            assert np.isfinite(kf.P).all(), "Covariance contains NaN or Inf after predict"
+
+            kf.update(np.array([[np.random.randn()]]))
+            assert np.isfinite(kf.P).all(), "Covariance contains NaN or Inf after update"
+
+    @pytest.mark.property
+    def test_state_shape_consistency(self, kf_2d_fixture):
+        """Property: State shape must remain constant throughout filter lifecycle."""
+        kf = kf_2d_fixture
+        initial_shape = kf.x.shape
+
+        for _ in range(20):
+            kf.predict()
+            assert kf.x.shape == initial_shape, "State shape changed after predict"
+
+            kf.update(np.array([[0.0], [0.0]]))
+            assert kf.x.shape == initial_shape, "State shape changed after update"
+
+    @pytest.mark.property
+    def test_covariance_shape_consistency(self, kf_2d_fixture):
+        """Property: Covariance shape must remain constant."""
+        kf = kf_2d_fixture
+        initial_shape = kf.P.shape
+
+        for _ in range(20):
+            kf.predict()
+            assert kf.P.shape == initial_shape, "Covariance shape changed"
+
+            kf.update(np.array([[0.0], [0.0]]))
+            assert kf.P.shape == initial_shape, "Covariance shape changed"
+
+    @pytest.mark.property
+    def test_kalman_gain_valid_dimensions(self, kf_1d_fixture):
+        """Property: Kalman gain K must have correct dimensions."""
+        kf = kf_1d_fixture
+
+        for _ in range(10):
+            kf.predict()
+            kf.update(np.array([[0.0]]))
+
+            # K should be dim_x × dim_z
+            assert kf.K.shape == (kf.dim_x, kf.dim_z)
+
+    @pytest.mark.property
+    @given(
+        scale_q=st.floats(min_value=0.001, max_value=10.0),
+        scale_r=st.floats(min_value=0.001, max_value=10.0),
+    )
+    def test_uncertainty_growth_with_process_noise(self, scale_q, scale_r):
+        """Property: Larger Q should result in larger steady-state covariance."""
+        dt = 0.1
+
+        def make_filter(q_scale, r_scale):
+            kf = KalmanFilter(dim_x=2, dim_z=1)
+            kf.x = np.array([[0.0], [1.0]])
+            kf.F = np.array([[1.0, dt], [0.0, 1.0]])
+            kf.H = np.array([[1.0, 0.0]])
+            kf.P = np.eye(2) * 100.0
+            kf.Q = Q_discrete_white_noise(dim=2, dt=dt, var=0.1 * q_scale)
+            kf.R = np.array([[1.0 * r_scale]])
+            return kf
+
+        # Run to steady state
+        kf_low_q = make_filter(0.1, scale_r)
+        kf_high_q = make_filter(1.0, scale_r)
+
+        for _ in range(50):
+            for kf in [kf_low_q, kf_high_q]:
+                kf.predict()
+                kf.update(np.array([[0.0]]))
+
+        # Higher Q should lead to higher uncertainty
+        trace_low_q = np.trace(kf_low_q.P)
+        trace_high_q = np.trace(kf_high_q.P)
+        assert trace_high_q > trace_low_q
+
+
+# ============================================================================
+# PART 3: EDGE CASE AND ERROR HANDLING TESTS
+# ============================================================================
+
+
+class TestKalmanFilterEdgeCases:
+    """Test edge cases and error handling."""
+
+    @pytest.mark.unit
+    def test_zero_measurement_noise(self):
+        """Test handling of zero measurement noise (R=0)."""
+        kf = KalmanFilter(dim_x=2, dim_z=1)
+
+        kf.x = np.array([[0.0], [1.0]])
+        kf.F = np.array([[1.0, 0.1], [0.0, 1.0]])
+        kf.H = np.array([[1.0, 0.0]])
+        kf.P = np.eye(2) * 100.0
+        kf.Q = np.zeros((2, 2))
+        kf.R = np.zeros((1, 1))
+
+        # Filter should not crash with R=0
+        kf.predict()
+        # Use a try-except to handle potential numerical issues
+        try:
+            kf.update(np.array([[1.0]]))
+            # If update succeeds, the state should be updated
+            assert np.isfinite(kf.x).all()
+        except np.linalg.LinAlgError:
+            # This is acceptable for zero measurement noise
+            pass
+
+    @pytest.mark.unit
+    def test_zero_process_noise(self):
+        """Test handling of zero process noise (Q=0)."""
+        kf = KalmanFilter(dim_x=2, dim_z=1)
+
+        kf.x = np.array([[0.0], [1.0]])
+        kf.F = np.array([[1.0, 0.1], [0.0, 1.0]])
+        kf.H = np.array([[1.0, 0.0]])
+        kf.P = np.eye(2) * 100.0
+        kf.Q = np.zeros((2, 2))
+        kf.R = np.array([[1.0]])
+
+        # Should handle Q=0 gracefully
+        for _ in range(5):
+            kf.predict()
+            kf.update(np.array([[0.0]]))
+
+            assert np.isfinite(kf.x).all()
+            assert np.isfinite(kf.P).all()
+
+    @pytest.mark.unit
+    def test_very_small_measurement_noise(self):
+        """Test handling of very small measurement noise."""
+        kf = KalmanFilter(dim_x=2, dim_z=1)
+
+        kf.x = np.array([[0.0], [1.0]])
+        kf.F = np.array([[1.0, 0.1], [0.0, 1.0]])
+        kf.H = np.array([[1.0, 0.0]])
+        kf.P = np.eye(2) * 100.0
+        kf.Q = Q_discrete_white_noise(dim=2, dt=0.1, var=0.01)
+        kf.R = np.array([[1e-10]])
+
+        for _ in range(10):
+            kf.predict()
+            kf.update(np.array([[0.0]]))
+
+            assert np.isfinite(kf.x).all()
+            assert np.isfinite(kf.P).all()
+            assert_matrix_psd(kf.P)
+
+    @pytest.mark.unit
+    def test_large_measurement_residual(self):
+        """Test handling of outlier measurements."""
+        kf = KalmanFilter(dim_x=2, dim_z=1)
+
+        kf.x = np.array([[0.0], [1.0]])
+        kf.F = np.array([[1.0, 0.1], [0.0, 1.0]])
+        kf.H = np.array([[1.0, 0.0]])
+        kf.P = np.eye(2) * 100.0
+        kf.Q = Q_discrete_white_noise(dim=2, dt=0.1, var=0.01)
+        kf.R = np.array([[1.0]])
+
+        # Update with a huge outlier
+        kf.predict()
+        kf.update(np.array([[1000.0]]))
+
+        # Filter should still be valid
+        assert np.isfinite(kf.x).all()
+        assert np.isfinite(kf.P).all()
+        assert_matrix_psd(kf.P)
+
+    @pytest.mark.unit
+    def test_repeated_predictions_without_update(self):
+        """Test multiple predictions without measurements."""
+        kf = KalmanFilter(dim_x=2, dim_z=1)
+
+        kf.x = np.array([[0.0], [1.0]])
+        kf.F = np.array([[1.0, 0.1], [0.0, 1.0]])
+        kf.H = np.array([[1.0, 0.0]])
+        kf.P = np.eye(2) * 10.0
+        kf.Q = Q_discrete_white_noise(dim=2, dt=0.1, var=0.01)
+        kf.R = np.array([[1.0]])
+
+        for _ in range(20):
+            kf.predict()
+            assert np.isfinite(kf.x).all()
+            assert np.isfinite(kf.P).all()
+
+    @pytest.mark.unit
+    def test_single_dimension_filter(self):
+        """Test degenerate case of 1x1 filter."""
+        kf = KalmanFilter(dim_x=1, dim_z=1)
+
+        kf.x = np.array([[0.0]])
+        kf.F = np.array([[1.0]])
+        kf.H = np.array([[1.0]])
+        kf.P = np.array([[100.0]])
+        kf.Q = np.array([[0.01]])
+        kf.R = np.array([[1.0]])
+
+        for _ in range(10):
+            kf.predict()
+            kf.update(np.array([[np.random.randn()]]))
+
+            assert kf.x.shape == (1, 1)
+            assert kf.P.shape == (1, 1)
+
+    @pytest.mark.unit
+    def test_high_dimensional_filter(self):
+        """Test high-dimensional filter (10 states)."""
+        dim = 10
+        kf = KalmanFilter(dim_x=dim, dim_z=dim)
+
+        kf.x = np.zeros((dim, 1))
+        kf.F = np.eye(dim)
+        kf.H = np.eye(dim)
+        kf.P = np.eye(dim) * 100.0
+        kf.Q = np.eye(dim) * 0.01
+        kf.R = np.eye(dim)
+
+        for _ in range(5):
+            kf.predict()
+            z = np.random.randn(dim, 1)
+            kf.update(z)
+
+            assert kf.x.shape == (dim, 1)
+            assert kf.P.shape == (dim, dim)
+
+    @pytest.mark.unit
+    def test_mismatched_measurement_dimensions(self):
+        """Test error handling for measurement dimension mismatch."""
+        kf = KalmanFilter(dim_x=2, dim_z=1)
+
+        kf.x = np.array([[0.0], [1.0]])
+        kf.F = np.array([[1.0, 0.1], [0.0, 1.0]])
+        kf.H = np.array([[1.0, 0.0]])
+        kf.P = np.eye(2) * 100.0
+        kf.Q = np.zeros((2, 2))
+        kf.R = np.array([[1.0]])
+
+        kf.predict()
+
+        # Provide wrong measurement dimension
+        with pytest.raises((ValueError, AssertionError)):
+            kf.update(np.array([[1.0], [2.0]]))  # Should be 1D, not 2D
+
+    @pytest.mark.unit
+    def test_very_large_state_values(self):
+        """Test handling of very large state values."""
+        kf = KalmanFilter(dim_x=2, dim_z=1)
+
+        kf.x = np.array([[1e6], [1.0]])
+        kf.F = np.array([[1.0, 0.1], [0.0, 1.0]])
+        kf.H = np.array([[1.0, 0.0]])
+        kf.P = np.eye(2) * 100.0
+        kf.Q = Q_discrete_white_noise(dim=2, dt=0.1, var=0.01)
+        kf.R = np.array([[1.0]])
+
+        # Filter should handle large values
+        for _ in range(10):
+            kf.predict()
+            kf.update(np.array([[1e6 + np.random.randn()]]))
+
+            assert np.isfinite(kf.x).all()
+            assert np.isfinite(kf.P).all()
+
+
+# ============================================================================
+# PART 4: NUMERICAL STABILITY TESTS
+# ============================================================================
+
+
+class TestKalmanFilterNumericalStability:
+    """Test numerical stability over long runs and ill-conditioned systems."""
+
+    @pytest.mark.slow
+    @pytest.mark.numerical
+    def test_long_running_stability_100_iterations(self):
+        """Test stability over 100 iterations."""
+        kf = KalmanFilter(dim_x=2, dim_z=1)
+
+        kf.x = np.array([[0.0], [1.0]])
+        kf.F = np.array([[1.0, 0.1], [0.0, 1.0]])
+        kf.H = np.array([[1.0, 0.0]])
+        kf.P = np.eye(2) * 100.0
+        kf.Q = Q_discrete_white_noise(dim=2, dt=0.1, var=0.1)
+        kf.R = np.array([[1.0]])
+
+        np.random.seed(42)
+
+        for _ in range(100):
+            kf.predict()
+            z = np.random.randn() * np.sqrt(1.0) + np.sin(_)
+            kf.update(np.array([[z]]))
+
+            # Check numerical health at each iteration
+            assert np.isfinite(kf.x).all(), f"NaN/Inf in x at iteration {_}"
+            assert np.isfinite(kf.P).all(), f"NaN/Inf in P at iteration {_}"
+            assert_matrix_psd(kf.P)
+
+    @pytest.mark.slow
+    @pytest.mark.numerical
+    def test_long_running_stability_1000_iterations(self):
+        """Test stability over 1000 iterations."""
+        kf = KalmanFilter(dim_x=2, dim_z=1)
+
+        kf.x = np.array([[0.0], [1.0]])
+        kf.F = np.array([[1.0, 0.1], [0.0, 1.0]])
+        kf.H = np.array([[1.0, 0.0]])
+        kf.P = np.eye(2) * 10.0
+        kf.Q = Q_discrete_white_noise(dim=2, dt=0.1, var=0.05)
+        kf.R = np.array([[1.0]])
+
+        np.random.seed(42)
+
+        traces = []
+        for i in range(1000):
+            kf.predict()
+            z = 0.5 * np.sin(i * 0.01) + np.random.randn() * 0.1
+            kf.update(np.array([[z]]))
+
+            traces.append(np.trace(kf.P))
+
+            if i % 100 == 0:
+                # Periodically verify health
+                assert np.isfinite(kf.P).all()
+                assert_matrix_psd(kf.P)
+
+        # Covariance should converge, not diverge
+        late_mean = np.mean(traces[-100:])
+        assert not np.isnan(late_mean)
+        assert not np.isinf(late_mean)
+
+    @pytest.mark.numerical
+    def test_ill_conditioned_system(self):
+        """Test filter stability with ill-conditioned system.
+
+        Setup: State variables with very different magnitudes.
+        Example: Position in meters (~1) and orientation in microradians (~1e-6).
+        """
+        kf = KalmanFilter(dim_x=2, dim_z=1)
+
+        # State: [position_meters, orientation_microradians]
+        kf.x = np.array([[100.0], [0.000001]])
+        kf.F = np.array([[1.0, 0.1], [0.0, 1.0]])
+        kf.H = np.array([[1.0, 0.0]])
+        kf.P = np.diag([100.0, 1e-12])  # Very different scales
+        kf.Q = np.diag([0.01, 1e-14])
+        kf.R = np.array([[1.0]])
+
+        for _ in range(20):
+            kf.predict()
+            kf.update(np.array([[100.0 + np.random.randn()]]))
+
+            # Should remain stable despite ill conditioning
+            assert np.isfinite(kf.x).all()
+            assert np.isfinite(kf.P).all()
+            assert_matrix_psd(kf.P)
+
+    @pytest.mark.numerical
+    def test_filter_with_changing_measurement_noise(self):
+        """Test adaptive filter with changing measurement noise."""
+        kf = KalmanFilter(dim_x=2, dim_z=1)
+
+        kf.x = np.array([[0.0], [1.0]])
+        kf.F = np.array([[1.0, 0.1], [0.0, 1.0]])
+        kf.H = np.array([[1.0, 0.0]])
+        kf.P = np.eye(2) * 100.0
+        kf.Q = Q_discrete_white_noise(dim=2, dt=0.1, var=0.01)
+
+        np.random.seed(42)
+
+        for i in range(50):
+            kf.predict()
+
+            # Measurement noise increases over time
+            r = 1.0 + 0.01 * i
+            kf.R = np.array([[r**2]])
+
+            z = i * 0.1 + np.random.randn() * r
+            kf.update(np.array([[z]]))
+
+            assert np.isfinite(kf.x).all()
+            assert np.isfinite(kf.P).all()
+
+    @pytest.mark.numerical
+    def test_filter_convergence_rate(self):
+        """Test that filter converges at expected rate."""
+        kf = KalmanFilter(dim_x=2, dim_z=1)
+
+        kf.x = np.array([[0.0], [1.0]])
+        kf.F = np.array([[1.0, 0.1], [0.0, 1.0]])
+        kf.H = np.array([[1.0, 0.0]])
+        kf.P = np.eye(2) * 1000.0
+        kf.Q = Q_discrete_white_noise(dim=2, dt=0.1, var=0.01)
+        kf.R = np.array([[1.0]])
+
+        traces = []
+        for _ in range(50):
+            kf.predict()
+            kf.update(np.array([[0.0]]))
+            traces.append(np.trace(kf.P))
+
+        # Filter should converge monotonically (with some numerical tolerance)
+        for i in range(1, len(traces)):
+            assert traces[i] <= traces[i - 1] + 1e-6
diff --git a/bayesian_filters/testing_utils/__init__.py b/bayesian_filters/testing_utils/__init__.py
new file mode 100644
index 0000000..f3a778e
--- /dev/null
+++ b/bayesian_filters/testing_utils/__init__.py
@@ -0,0 +1,42 @@
+"""Testing utilities for bayesian_filters.
+
+This package provides fixtures, helpers, and analytical solutions for testing
+Kalman filters and related algorithms.
+
+Submodules:
+- fixtures: Common pytest fixtures for filter testing
+- numerical: Analytical solutions and reference implementations
+- helpers: Utility functions for test development
+"""
+
+from bayesian_filters.testing_utils.fixtures import (
+    kf_1d,
+    kf_2d,
+    ekf_2d,
+    ukf_3d,
+    sensor_sim,
+    noisy_measurements,
+)
+from bayesian_filters.testing_utils.numerical import (
+    ConstantVelocitySolution,
+    ConstantAccelerationSolution,
+)
+from bayesian_filters.testing_utils.helpers import (
+    assert_matrix_psd,
+    assert_matrix_symmetric,
+    assert_filter_stable,
+)
+
+__all__ = [
+    "kf_1d",
+    "kf_2d",
+    "ekf_2d",
+    "ukf_3d",
+    "sensor_sim",
+    "noisy_measurements",
+    "ConstantVelocitySolution",
+    "ConstantAccelerationSolution",
+    "assert_matrix_psd",
+    "assert_matrix_symmetric",
+    "assert_filter_stable",
+]
diff --git a/bayesian_filters/testing_utils/fixtures.py b/bayesian_filters/testing_utils/fixtures.py
new file mode 100644
index 0000000..620c21b
--- /dev/null
+++ b/bayesian_filters/testing_utils/fixtures.py
@@ -0,0 +1,234 @@
+"""Common pytest fixtures for Bayesian filter testing."""
+
+import numpy as np
+from bayesian_filters.kalman import (
+    KalmanFilter,
+    ExtendedKalmanFilter,
+    UnscentedKalmanFilter,
+    MerweScaledSigmaPoints,
+)
+from bayesian_filters.common import Q_discrete_white_noise
+
+
+class SimpleSensorSim:
+    """Simple 1D sensor simulator."""
+
+    def __init__(self, pos=0.0, vel=1.0, dt=0.1, noise_std=1.0):
+        """Initialize sensor simulator.
+
+        Parameters
+        ----------
+        pos : float
+            Initial position
+        vel : float
+            Velocity (units per time step)
+        dt : float
+            Time step
+        noise_std : float
+            Measurement noise standard deviation
+        """
+        self.pos = pos
+        self.vel = vel
+        self.dt = dt
+        self.noise_std = noise_std
+
+    def read(self):
+        """Read next measurement with noise."""
+        self.pos += self.vel * self.dt
+        return self.pos + np.random.normal(0, self.noise_std)
+
+    def true_position(self):
+        """Return current true position."""
+        return self.pos
+
+
+def kf_1d():
+    """Create a 1D constant-velocity Kalman filter fixture.
+
+    Returns
+    -------
+    KalmanFilter
+        1D Kalman filter with state [position, velocity]
+    """
+    dt = 0.1
+    kf = KalmanFilter(dim_x=2, dim_z=1)
+
+    # State transition matrix (constant velocity model)
+    kf.F = np.array([[1.0, dt], [0.0, 1.0]])
+
+    # Measurement matrix (measure position only)
+    kf.H = np.array([[1.0, 0.0]])
+
+    # Initial covariance
+    kf.P = np.eye(2) * 100.0
+
+    # Measurement noise
+    kf.R = np.array([[1.0]])
+
+    # Process noise
+    kf.Q = Q_discrete_white_noise(dim=2, dt=dt, var=0.1)
+
+    # Initial state
+    kf.x = np.array([[0.0], [1.0]])
+
+    return kf
+
+
+def kf_2d():
+    """Create a 2D constant-velocity Kalman filter fixture.
+
+    Returns
+    -------
+    KalmanFilter
+        2D Kalman filter with state [x, vx, y, vy]
+    """
+    dt = 0.1
+    kf = KalmanFilter(dim_x=4, dim_z=2)
+
+    # State transition matrix
+    F = np.eye(4)
+    F[0, 1] = dt
+    F[2, 3] = dt
+    kf.F = F
+
+    # Measurement matrix (measure both positions)
+    kf.H = np.array([[1.0, 0.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0]])
+
+    # Initial covariance
+    kf.P = np.eye(4) * 100.0
+
+    # Measurement noise
+    kf.R = np.eye(2)
+
+    # Process noise
+    q = Q_discrete_white_noise(dim=2, dt=dt, var=0.1)
+    kf.Q = np.zeros((4, 4))
+    kf.Q[:2, :2] = q
+    kf.Q[2:, 2:] = q
+
+    # Initial state
+    kf.x = np.array([[0.0], [0.0], [0.0], [0.0]])
+
+    return kf
+
+
+def ekf_2d():
+    """Create a 2D Extended Kalman filter fixture.
+
+    Returns
+    -------
+    ExtendedKalmanFilter
+        2D EKF for nonlinear measurement function
+    """
+    dt = 0.1
+    ekf = ExtendedKalmanFilter(dim_x=2, dim_z=1)
+
+    # State transition matrix
+    ekf.F = np.array([[1.0, dt], [0.0, 1.0]])
+
+    # State transition function
+    def fx(x, dt):
+        return np.dot(ekf.F, x)
+
+    # Measurement function (range to origin)
+    def hx(x):
+        return np.sqrt(x[0] ** 2 + x[1] ** 2)
+
+    # Jacobian of measurement function
+    def H(x):
+        r = np.sqrt(x[0] ** 2 + x[1] ** 2)
+        if r < 1e-10:
+            r = 1e-10
+        return np.array([[x[0] / r, x[1] / r]])
+
+    ekf.fx = fx
+    ekf.hx = hx
+    ekf.H = H
+
+    # Covariances
+    ekf.P = np.eye(2) * 100.0
+    ekf.R = np.array([[1.0]])
+    ekf.Q = Q_discrete_white_noise(dim=2, dt=dt, var=0.1)
+
+    # Initial state
+    ekf.x = np.array([10.0, 0.0])
+
+    return ekf
+
+
+def ukf_3d():
+    """Create a 3D Unscented Kalman filter fixture.
+
+    Returns
+    -------
+    UnscentedKalmanFilter
+        3D UKF with standard configuration
+    """
+    dt = 0.1
+
+    # Create sigma points
+    sigmas = MerweScaledSigmaPoints(n=3, alpha=0.1, beta=2.0, kappa=0.0)
+
+    ukf = UnscentedKalmanFilter(
+        dim_x=3,
+        dim_z=3,
+        dt=dt,
+        hx=lambda x: x,  # Identity measurement function
+        fx=lambda x, dt: x,  # Identity process function
+        points=sigmas,
+    )
+
+    # Covariances
+    ukf.P = np.eye(3) * 100.0
+    ukf.R = np.eye(3)
+    ukf.Q = Q_discrete_white_noise(dim=3, dt=dt, var=0.1)
+
+    # Initial state
+    ukf.x = np.array([0.0, 0.0, 0.0])
+
+    return ukf
+
+
+def sensor_sim(pos=0.0, vel=1.0, dt=0.1, noise_std=1.0):
+    """Create a simple sensor simulator fixture.
+
+    Parameters
+    ----------
+    pos : float
+        Initial position
+    vel : float
+        Velocity
+    dt : float
+        Time step
+    noise_std : float
+        Measurement noise standard deviation
+
+    Returns
+    -------
+    SimpleSensorSim
+        Configured sensor simulator
+    """
+    return SimpleSensorSim(pos=pos, vel=vel, dt=dt, noise_std=noise_std)
+
+
+def noisy_measurements(ground_truth, noise_std, seed=42):
+    """Generate noisy measurements from ground truth.
+
+    Parameters
+    ----------
+    ground_truth : array_like
+        True values
+    noise_std : float
+        Standard deviation of measurement noise
+    seed : int
+        Random seed for reproducibility
+
+    Returns
+    -------
+    np.ndarray
+        Measurements with added Gaussian noise
+    """
+    np.random.seed(seed)
+    ground_truth = np.asarray(ground_truth)
+    noise = np.random.normal(0, noise_std, size=ground_truth.shape)
+    return ground_truth + noise
diff --git a/bayesian_filters/testing_utils/helpers.py b/bayesian_filters/testing_utils/helpers.py
new file mode 100644
index 0000000..be23e79
--- /dev/null
+++ b/bayesian_filters/testing_utils/helpers.py
@@ -0,0 +1,218 @@
+"""Helper functions for testing Bayesian filters."""
+
+import numpy as np
+
+
+def assert_matrix_psd(M, rtol=1e-5, atol=1e-8, name="Matrix"):
+    """Assert that matrix is positive semi-definite.
+
+    Parameters
+    ----------
+    M : np.ndarray
+        Matrix to test
+    rtol : float
+        Relative tolerance for eigenvalue checks
+    atol : float
+        Absolute tolerance for eigenvalue checks
+    name : str
+        Name of matrix for error messages
+
+    Raises
+    ------
+    AssertionError
+        If matrix is not positive semi-definite
+    """
+    M = np.asarray(M)
+
+    # Check symmetry
+    if not np.allclose(M, M.T, rtol=rtol, atol=atol):
+        raise AssertionError(f"{name} is not symmetric")
+
+    # Check eigenvalues
+    eigvals = np.linalg.eigvalsh(M)
+    min_eigval = np.min(eigvals)
+
+    if min_eigval < -atol:
+        raise AssertionError(f"{name} is not positive semi-definite. Minimum eigenvalue: {min_eigval}")
+
+
+def assert_matrix_symmetric(M, rtol=1e-5, atol=1e-8, name="Matrix"):
+    """Assert that matrix is symmetric.
+
+    Parameters
+    ----------
+    M : np.ndarray
+        Matrix to test
+    rtol : float
+        Relative tolerance
+    atol : float
+        Absolute tolerance
+    name : str
+        Name of matrix for error messages
+
+    Raises
+    ------
+    AssertionError
+        If matrix is not symmetric
+    """
+    M = np.asarray(M)
+
+    if not np.allclose(M, M.T, rtol=rtol, atol=atol):
+        diff = np.max(np.abs(M - M.T))
+        raise AssertionError(f"{name} is not symmetric. Max difference: {diff}")
+
+
+def assert_filter_stable(P, max_variance=1e10):
+    """Assert that filter covariance matrix is stable.
+
+    Checks that diagonal elements don't explode.
+
+    Parameters
+    ----------
+    P : np.ndarray
+        Filter covariance matrix
+    max_variance : float
+        Maximum allowed variance on diagonal
+
+    Raises
+    ------
+    AssertionError
+        If filter appears unstable
+    """
+    P = np.asarray(P)
+    diag = np.diag(P)
+
+    if np.any(diag > max_variance):
+        raise AssertionError(f"Filter covariance is unstable. Max diagonal: {np.max(diag)}")
+
+    if np.any(np.isnan(diag)):
+        raise AssertionError("Filter covariance contains NaN values")
+
+    if np.any(np.isinf(diag)):
+        raise AssertionError("Filter covariance contains infinite values")
+
+
+def compute_condition_number(M):
+    """Compute condition number of matrix.
+
+    Parameters
+    ----------
+    M : np.ndarray
+        Square matrix
+
+    Returns
+    -------
+    float
+        Condition number (ratio of largest to smallest singular value)
+    """
+    M = np.asarray(M)
+    s = np.linalg.svd(M, compute_uv=False)
+    return np.max(s) / np.max(np.min(s), 1e-15)
+
+
+def assert_well_conditioned(M, max_condition=1e6, name="Matrix"):
+    """Assert that matrix is well-conditioned.
+
+    Parameters
+    ----------
+    M : np.ndarray
+        Matrix to test
+    max_condition : float
+        Maximum allowed condition number
+    name : str
+        Name of matrix for error messages
+
+    Raises
+    ------
+    AssertionError
+        If matrix is ill-conditioned
+    """
+    cond = compute_condition_number(M)
+
+    if cond > max_condition:
+        raise AssertionError(f"{name} is ill-conditioned. Condition number: {cond:.2e}")
+
+
+def assert_convergence(errors, tolerance=0.1, rtol=0.05):
+    """Assert that error sequence converges.
+
+    Parameters
+    ----------
+    errors : array_like
+        Sequence of errors
+    tolerance : float
+        Final error should be below this
+    rtol : float
+        Relative tolerance for "no divergence"
+
+    Raises
+    ------
+    AssertionError
+        If sequence doesn't converge
+    """
+    errors = np.asarray(errors)
+
+    # Check final error
+    if errors[-1] > tolerance:
+        raise AssertionError(f"Errors did not converge. Final error: {errors[-1]}")
+
+    # Check for divergence in later part
+    final_half = errors[len(errors) // 2 :]
+    if np.max(final_half) > np.min(final_half) * (1 + rtol):
+        raise AssertionError(f"Errors appear to diverge. Min: {np.min(final_half)}, Max: {np.max(final_half)}")
+
+
+def mahalanobis_distance(x, mean, cov):
+    """Compute Mahalanobis distance.
+
+    Parameters
+    ----------
+    x : array_like
+        Vector
+    mean : array_like
+        Mean vector
+    cov : array_like
+        Covariance matrix
+
+    Returns
+    -------
+    float
+        Mahalanobis distance
+    """
+    x = np.asarray(x).flatten()
+    mean = np.asarray(mean).flatten()
+    cov = np.asarray(cov)
+
+    diff = x - mean
+    try:
+        cov_inv = np.linalg.inv(cov)
+    except np.linalg.LinAlgError:
+        cov_inv = np.linalg.pinv(cov)
+
+    return np.sqrt(np.dot(diff, np.dot(cov_inv, diff)))
+
+
+def normalized_error(estimate, truth):
+    """Compute normalized estimation error.
+
+    Parameters
+    ----------
+    estimate : array_like
+        Estimated values
+    truth : array_like
+        True values
+
+    Returns
+    -------
+    np.ndarray
+        Normalized errors
+    """
+    estimate = np.asarray(estimate).flatten()
+    truth = np.asarray(truth).flatten()
+
+    # Avoid division by zero
+    with np.errstate(divide="ignore", invalid="ignore"):
+        error = (estimate - truth) / np.abs(truth)
+        error = np.nan_to_num(error, nan=0.0, posinf=0.0, neginf=0.0)
+
+    return error
diff --git a/bayesian_filters/testing_utils/numerical.py b/bayesian_filters/testing_utils/numerical.py
new file mode 100644
index 0000000..80fbdda
--- /dev/null
+++ b/bayesian_filters/testing_utils/numerical.py
@@ -0,0 +1,230 @@
+"""Analytical solutions and reference implementations for testing."""
+
+import numpy as np
+
+
+class ConstantVelocitySolution:
+    """Analytical solution for constant velocity filter.
+
+    For a 1D system with constant velocity:
+    - State: [position, velocity]
+    - Process model: x(k+1) = x(k) + v(k)*dt
+    - Measurement: position only
+
+    With zero process noise and zero initial velocity error,
+    the filter should track perfectly.
+    """
+
+    @staticmethod
+    def position_at_time(t, x0=0.0, v0=1.0):
+        """Analytical position at time t.
+
+        Parameters
+        ----------
+        t : float
+            Time
+        x0 : float
+            Initial position
+        v0 : float
+            Constant velocity
+
+        Returns
+        -------
+        float
+            Position at time t
+        """
+        return x0 + v0 * t
+
+    @staticmethod
+    def steady_state_error(R, Q, dt, var_init=None):
+        """Calculate steady-state estimation error.
+
+        Parameters
+        ----------
+        R : float
+            Measurement noise variance
+        Q : float
+            Process noise variance per dt
+        dt : float
+            Time step
+        var_init : float, optional
+            Initial position variance
+
+        Returns
+        -------
+        float
+            Steady-state position error variance
+        """
+        # Simplified steady-state for constant velocity model
+        # With proper tuning, error should converge to sqrt(Q*R)
+        return np.sqrt(Q * R / dt)
+
+
+class ConstantAccelerationSolution:
+    """Analytical solution for constant acceleration filter.
+
+    For a 1D system with constant acceleration:
+    - State: [position, velocity, acceleration]
+    - Constant acceleration assumption
+    """
+
+    @staticmethod
+    def position_at_time(t, x0=0.0, v0=0.0, a=1.0):
+        """Analytical position at time t.
+
+        Parameters
+        ----------
+        t : float
+            Time
+        x0 : float
+            Initial position
+        v0 : float
+            Initial velocity
+        a : float
+            Constant acceleration
+
+        Returns
+        -------
+        float
+            Position at time t
+        """
+        return x0 + v0 * t + 0.5 * a * t**2
+
+    @staticmethod
+    def velocity_at_time(t, v0=0.0, a=1.0):
+        """Analytical velocity at time t.
+
+        Parameters
+        ----------
+        t : float
+            Time
+        v0 : float
+            Initial velocity
+        a : float
+            Constant acceleration
+
+        Returns
+        -------
+        float
+            Velocity at time t
+        """
+        return v0 + a * t
+
+
+class LinearSystemSolution:
+    """Analytical solution for linear systems.
+
+    Computes the optimal state estimate for a linear system with
+    known process and measurement models.
+    """
+
+    @staticmethod
+    def kalman_gain(P, H, R):
+        """Compute Kalman gain.
+
+        Parameters
+        ----------
+        P : np.ndarray
+            Prediction covariance (n x n)
+        H : np.ndarray
+            Measurement matrix (m x n)
+        R : np.ndarray
+            Measurement noise covariance (m x m)
+
+        Returns
+        -------
+        np.ndarray
+            Kalman gain (n x m)
+        """
+        S = np.dot(H, np.dot(P, H.T)) + R  # Innovation covariance
+        K = np.dot(P, np.dot(H.T, np.linalg.inv(S)))
+        return K
+
+    @staticmethod
+    def innovation_covariance(P, H, R):
+        """Compute innovation (measurement residual) covariance.
+
+        Parameters
+        ----------
+        P : np.ndarray
+            Prediction covariance (n x n)
+        H : np.ndarray
+            Measurement matrix (m x n)
+        R : np.ndarray
+            Measurement noise covariance (m x m)
+
+        Returns
+        -------
+        np.ndarray
+            Innovation covariance (m x m)
+        """
+        return np.dot(H, np.dot(P, H.T)) + R
+
+    @staticmethod
+    def updated_covariance(K, H, P):
+        """Compute updated state covariance after Kalman update.
+
+        Parameters
+        ----------
+        K : np.ndarray
+            Kalman gain (n x m)
+        H : np.ndarray
+            Measurement matrix (m x n)
+        P : np.ndarray
+            Prediction covariance (n x n)
+
+        Returns
+        -------
+        np.ndarray
+            Updated covariance (n x n)
+        """
+        return (np.eye(P.shape[0]) - np.dot(K, H)) @ P
+
+
+class WhiteNoiseAcceleration:
+    """White noise acceleration process model.
+
+    Reference: Bar-Shalom et al., "Estimation with Applications
+    to Tracking and Navigation"
+    """
+
+    @staticmethod
+    def discrete_process_noise(var, dt, order=2):
+        """Compute discrete process noise covariance.
+
+        For continuous white noise acceleration model with
+        discrete sampling at interval dt.
+
+        Parameters
+        ----------
+        var : float
+            Acceleration variance (continuous)
+        dt : float
+            Sampling interval
+        order : int, default=2
+            Model order (2=const velocity, 3=const acceleration)
+
+        Returns
+        -------
+        np.ndarray
+            Process noise covariance matrix
+        """
+        if order == 2:
+            # Constant velocity model
+            return var * np.array(
+                [
+                    [dt**3 / 3, dt**2 / 2],
+                    [dt**2 / 2, dt],
+                ]
+            )
+        elif order == 3:
+            # Constant acceleration model
+            return var * np.array(
+                [
+                    [dt**5 / 20, dt**4 / 8, dt**3 / 6],
+                    [dt**4 / 8, dt**3 / 3, dt**2 / 2],
+                    [dt**3 / 6, dt**2 / 2, dt],
+                ]
+            )
+        else:
+            raise ValueError(f"Unsupported order: {order}")
diff --git a/conftest.py b/conftest.py
new file mode 100644
index 0000000..249d74e
--- /dev/null
+++ b/conftest.py
@@ -0,0 +1,73 @@
+"""Root pytest configuration for bayesian_filters."""
+
+import pytest
+import numpy as np
+from bayesian_filters.testing_utils.fixtures import (
+    kf_1d,
+    kf_2d,
+    ekf_2d,
+    ukf_3d,
+    sensor_sim,
+    noisy_measurements,
+)
+
+
+# Export fixtures so they're available to all tests
+pytest_plugins = ["pytest_benchmark"]
+
+
+@pytest.fixture
+def kf_1d_fixture():
+    """1D Kalman filter fixture."""
+    return kf_1d()
+
+
+@pytest.fixture
+def kf_2d_fixture():
+    """2D Kalman filter fixture."""
+    return kf_2d()
+
+
+@pytest.fixture
+def ekf_2d_fixture():
+    """2D Extended Kalman filter fixture."""
+    return ekf_2d()
+
+
+@pytest.fixture
+def ukf_3d_fixture():
+    """3D Unscented Kalman filter fixture."""
+    return ukf_3d()
+
+
+@pytest.fixture
+def sensor_sim_fixture():
+    """Sensor simulator fixture."""
+    return sensor_sim()
+
+
+def pytest_configure(config):
+    """Configure pytest."""
+    # Register custom markers
+    config.addinivalue_line("markers", "unit: mark test as a unit test")
+    config.addinivalue_line("markers", "integration: mark test as an integration test")
+    config.addinivalue_line("markers", "slow: mark test as slow (>1 second)")
+    config.addinivalue_line("markers", "benchmark: mark test as a benchmark")
+    config.addinivalue_line("markers", "property: mark test as property-based (hypothesis)")
+    config.addinivalue_line("markers", "numerical: mark test as numerical accuracy test")
+
+
+def pytest_collection_modifyitems(config, items):
+    """Modify collected items to add markers."""
+    for item in items:
+        # Auto-mark slow tests
+        if "benchmark" in item.nodeid:
+            item.add_marker(pytest.mark.benchmark)
+
+        # Auto-mark property tests
+        if "property" in item.nodeid:
+            item.add_marker(pytest.mark.property)
+
+
+# Set random seed for reproducibility
+np.random.seed(42)
diff --git a/pyproject.toml b/pyproject.toml
index 2e3090d..f3ab315 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -49,6 +49,10 @@ Documentation = "https://georgepearse.github.io/bayesian_filters"
 dev = [
     "pytest>=8.4.0",
     "pytest-cov>=7.0.0",
+    "hypothesis>=6.100.0",
+    "pytest-benchmark>=4.0.0",
+    "pytest-xdist>=3.5.0",
+    "pytest-mpl>=0.17.0",
     "zuban>=0.1.0",
 ]
 docs = [
@@ -68,7 +72,44 @@ testpaths = ["bayesian_filters"]
 python_files = ["test_*.py"]
 python_classes = ["Test*"]
 python_functions = ["test_*"]
-addopts = "-v --strict-markers"
+addopts = "-v --strict-markers -ra --tb=short"
+markers = [
+    "unit: Fast, isolated unit tests",
+    "integration: Tests multiple components together",
+    "slow: Tests that take >1 second to run",
+    "benchmark: Performance benchmark tests",
+    "property: Property-based tests (hypothesis)",
+    "numerical: Numerical accuracy tests",
+]
+norecursedirs = [".git", ".tox", "dist", "build", "*.egg", ".venv"]
+filterwarnings = [
+    "ignore::DeprecationWarning",
+    "ignore::PendingDeprecationWarning",
+]
+
+[tool.coverage.run]
+source = ["bayesian_filters"]
+relative_files = true
+omit = [
+    "*/__init__.py",
+    "*/tests/*",
+    "**/examples/**",
+]
+
+[tool.coverage.report]
+exclude_lines = [
+    "pragma: no cover",
+    "def __repr__",
+    "raise AssertionError",
+    "raise NotImplementedError",
+    "if __name__ == .__main__.:",
+    "if TYPE_CHECKING:",
+    "class .*\\bProtocol\\):",
+    "@(abc\\.)?abstractmethod",
+]
+skip_covered = false
+skip_empty = true
+precision = 2
 
 [tool.ruff]
 # Exclude test files and examples from most checks