.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_tutorials/06_uncertainty_estimation.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_tutorials_06_uncertainty_estimation.py: 06. Uncertainty Estimation ========================== This tutorial mainly explains the ``UncertaintyEstimator`` class. It demonstrates how ``UncertaintyEstimator`` converts source reconstruction outputs into calibration-ready uncertainty representations. It covers: - fixed-orientation aggregated marginal intervals; - free-orientation EEG aggregated ``marginal`` intervals; - free-orientation EEG aggregated ``full_cov`` ellipsoids; - pre-calibration empirical coverage curves before isotonic recalibration. .. GENERATED FROM PYTHON SOURCE LINES 20-36 Scientific motivation --------------------- Source estimation returns posterior means and posterior covariance matrices, but calibration does not operate on those objects directly. ``UncertaintyEstimator`` turns them into uncertainty representations that can be evaluated against ground truth. In the current workflow this means: - fixed orientation uses scalar ``posterior_var`` derived from the diagonal of the covariance; - free-orientation EEG can be evaluated either with pooled component-wise ``marginal`` intervals or with local 3D ``full_cov`` ellipsoids; - the default workflow uses temporally **aggregated** calibration, so means are averaged over time and covariance is scaled by ``1 / T``. .. GENERATED FROM PYTHON SOURCE LINES 36-52 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from mne.io.constants import FIFF from calibrain import ( SensorSimulator, SourceEstimator, SourceSimulator, UncertaintyEstimator, gamma_map_sflex, ) RANDOM_SEED = 53 .. GENERATED FROM PYTHON SOURCE LINES 53-66 Build a lightweight posterior example ------------------------------------- The tutorial is self-contained: simulate source activity, project it to EEG sensors, add noise, reconstruct sources with ``gamma_map_sflex``, then pass the posterior outputs into ``UncertaintyEstimator``. Units: - source amplitudes are in ``nAm``; - source coordinates for sFLEX are in ``m``; - EEG leadfields are interpreted as ``µV / nAm``; - sensor signals are therefore in ``µV``. .. GENERATED FROM PYTHON SOURCE LINES 66-115 .. code-block:: Python erp_config = { "tmin": -0.1, "tmax": 0.8, "stim_onset": 0.0, "sfreq": 100, "fmin": 2, "fmax": 8, "amplitude_distribution": { "median": 8.0, "sigma": 0.15, "clip": [2.0, 20.0], }, "random_erp_timing": False, "erp_min_length": 20, } nominal_coverages = np.linspace(0.0, 1.0, 11) uncertainty_estimator = UncertaintyEstimator(nominal_coverages=nominal_coverages) source_simulator = SourceSimulator(ERP_config=erp_config) sensor_simulator = SensorSimulator() times = np.arange(erp_config["tmin"], erp_config["tmax"], 1.0 / erp_config["sfreq"]) rng = np.random.default_rng(RANDOM_SEED) n_sensors = 16 n_sources = 32 src_coords = rng.normal(scale=0.04, size=(n_sources, 3)) leadfield_fixed = rng.normal(scale=0.03, size=(n_sensors, n_sources)) leadfield_fixed /= np.maximum( np.linalg.norm(leadfield_fixed, axis=0, keepdims=True), np.finfo(float).eps, ) leadfield_fixed *= 0.6 leadfield_free_eeg = rng.normal(scale=0.015, size=(n_sensors, n_sources, 3)) leadfield_free_eeg /= np.maximum( np.linalg.norm(leadfield_free_eeg, axis=0, keepdims=True), np.finfo(float).eps, ) leadfield_free_eeg *= 0.4 sensor_simulator.set_sensor_metadata( kind=FIFF.FIFFV_EEG_CH, units=FIFF.FIFF_UNIT_V, unitmult=FIFF.FIFF_UNITM_MU, coil_type=FIFF.FIFFV_COIL_EEG, ) .. GENERATED FROM PYTHON SOURCE LINES 116-121 Fixed orientation: from posterior covariance to ``posterior_var`` ----------------------------------------------------------------- For fixed orientation, calibration uses a 1D variance per source. This is the reduced representation written downstream into aggregated ``.npz`` files. .. GENERATED FROM PYTHON SOURCE LINES 121-155 .. code-block:: Python x_true_fixed, active_fixed = source_simulator.simulate( n_sources=n_sources, nnz=4, orientation_type="fixed", seed=RANDOM_SEED, ) y_fixed_clean, y_fixed_noisy, fixed_noise, fixed_eta = sensor_simulator.simulate( x=x_true_fixed, L=leadfield_fixed, alpha_SNR=0.7, sensor_white_noise_std=0.2, seed=RANDOM_SEED, ) fixed_noise_var = float(np.var(fixed_noise)) fixed_estimator = SourceEstimator( solver=gamma_map_sflex, solver_params={"max_iter": 150, "tol": 1e-7, "sigma": 0.01, "src_coords": src_coords}, noise_var=fixed_noise_var, n_orient=1, ) fixed_estimator.fit(leadfield_fixed, y_fixed_noisy) fixed_result = fixed_estimator.predict() posterior_var_fixed = uncertainty_estimator.posterior_variance_from_cov( fixed_result["posterior_cov"] ) print("fixed posterior_mean shape:", fixed_result["posterior_mean"].shape) print("fixed posterior_cov shape:", fixed_result["posterior_cov"].shape) print("fixed posterior_var shape:", posterior_var_fixed.shape) .. rst-class:: sphx-glr-script-out .. code-block:: none fixed posterior_mean shape: (32, 90) fixed posterior_cov shape: (32, 32) fixed posterior_var shape: (32,) .. GENERATED FROM PYTHON SOURCE LINES 156-161 Fixed orientation: aggregated intervals and pre-calibration curve ----------------------------------------------------------------- The active workflow uses aggregated calibration. ``UncertaintyEstimator`` averages source time courses over time and scales variance by ``1 / T``. .. GENERATED FROM PYTHON SOURCE LINES 161-177 .. code-block:: Python fixed_membership = uncertainty_estimator.aggregated_interval_membership( x_true=x_true_fixed, x_hat=fixed_result["posterior_mean"], posterior_var=posterior_var_fixed, nominal_coverage=0.9, ) fixed_curve = uncertainty_estimator.calibration_curve_intervals_aggregated( x_true=x_true_fixed, x_hat=fixed_result["posterior_mean"], posterior_var=posterior_var_fixed, ) print("fixed aggregated empirical coverage at 0.9:", fixed_membership["empirical_coverage"]) print("fixed interval_type:", fixed_curve["interval_type"]) .. rst-class:: sphx-glr-script-out .. code-block:: none fixed aggregated empirical coverage at 0.9: 1.0 fixed interval_type: marginal .. GENERATED FROM PYTHON SOURCE LINES 178-182 Plot a fixed-orientation aggregated interval example ---------------------------------------------------- The interval is shown for the time-averaged prediction of one active source. .. GENERATED FROM PYTHON SOURCE LINES 182-219 .. code-block:: Python fixed_source_idx = int(np.atleast_1d(active_fixed)[0]) fig, axes = plt.subplots(1, 2, figsize=(10, 4.2)) axes[0].axhline(fixed_membership["x_true_agg"][fixed_source_idx], color="black", label="true aggregated") axes[0].axhline(fixed_membership["x_hat_agg"][fixed_source_idx], color="C0", label="predicted aggregated") axes[0].fill_between( [0, 1], [fixed_membership["ci_lower"][fixed_source_idx]] * 2, [fixed_membership["ci_upper"][fixed_source_idx]] * 2, alpha=0.3, color="C0", label="90% interval", ) axes[0].set( xlim=(0, 1), xticks=[], ylabel="Aggregated source amplitude (nAm)", title=f"Fixed orientation: source {fixed_source_idx}", ) axes[0].legend(loc="best") axes[1].plot([0, 1], [0, 1], "--", color="0.5", label="perfect calibration") axes[1].plot( fixed_curve["nominal_coverages"], fixed_curve["empirical_coverages"], marker="o", label="fixed marginal", ) axes[1].set( xlabel="Nominal coverage", ylabel="Empirical coverage", title="Fixed orientation pre-calibration curve", ) axes[1].legend(loc="best") fig.tight_layout() .. image-sg:: /auto_tutorials/images/sphx_glr_06_uncertainty_estimation_001.png :alt: Fixed orientation: source 30, Fixed orientation pre-calibration curve :srcset: /auto_tutorials/images/sphx_glr_06_uncertainty_estimation_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 220-232 Free EEG orientation: ``marginal`` versus ``full_cov`` ------------------------------------------------------ For free-orientation EEG, calibration can use two uncertainty types: - ``marginal``: use only component-wise variances and pool over the three local orientation components; - ``full_cov``: use each local ``3 x 3`` covariance block and test coverage with 3D ellipsoids. Both use the same posterior mean and posterior covariance, but they answer slightly different questions. .. GENERATED FROM PYTHON SOURCE LINES 232-263 .. code-block:: Python x_true_free, active_free = source_simulator.simulate( n_sources=n_sources, nnz=4, orientation_type="free", coil_type=FIFF.FIFFV_COIL_EEG, seed=RANDOM_SEED + 1, ) y_free_clean, y_free_noisy, free_noise, free_eta = sensor_simulator.simulate( x=x_true_free, L=leadfield_free_eeg, alpha_SNR=0.7, sensor_white_noise_std=0.05, seed=RANDOM_SEED + 1, ) free_noise_var = float(np.var(free_noise)) free_estimator = SourceEstimator( solver=gamma_map_sflex, solver_params={"max_iter": 150, "tol": 1e-7, "sigma": 0.01, "src_coords": src_coords}, noise_var=free_noise_var, n_orient=3, ) free_estimator.fit(leadfield_free_eeg, y_free_noisy) free_result = free_estimator.predict() print("free EEG posterior_mean shape:", free_result["posterior_mean"].shape) print("free EEG posterior_mean_reshaped shape:", free_result["posterior_mean_reshaped"].shape) print("free EEG posterior_cov shape:", free_result["posterior_cov"].shape) .. rst-class:: sphx-glr-script-out .. code-block:: none free EEG posterior_mean shape: (96, 90) free EEG posterior_mean_reshaped shape: (32, 3, 90) free EEG posterior_cov shape: (96, 96) .. GENERATED FROM PYTHON SOURCE LINES 264-269 Compute aggregated pre-calibration curves ------------------------------------------------------------ ``marginal`` works with the same full covariance input, but only uses its diagonal entries source-by-source. ``full_cov`` uses the full local 3D blocks. .. GENERATED FROM PYTHON SOURCE LINES 269-300 .. code-block:: Python free_curve_marginal = uncertainty_estimator.calibration_curve_componentwise_eeg_free_aggregated( x_true=x_true_free, x_hat=free_result["posterior_mean_reshaped"], posterior_uncert=free_result["posterior_cov"], ) free_curve_full_cov = uncertainty_estimator.calibration_curve_ellipsoid_eeg_free_aggregated( x_true=x_true_free, x_hat=free_result["posterior_mean_reshaped"], posterior_cov=free_result["posterior_cov"], ) free_membership_marginal = uncertainty_estimator.aggregated_componentwise_interval_membership_free( x_true=x_true_free, x_hat=free_result["posterior_mean_reshaped"], posterior_uncert=free_result["posterior_cov"], nominal_coverage=0.9, n_orient=3, ) free_membership_full_cov = uncertainty_estimator.aggregated_ellipsoid_membership_eeg_free( x_true=x_true_free, x_hat=free_result["posterior_mean_reshaped"], posterior_cov=free_result["posterior_cov"], nominal_coverage=0.9, ) print("free marginal interval_type:", free_curve_marginal["interval_type"]) print("free full_cov interval_type:", free_curve_full_cov["interval_type"]) print("free marginal empirical coverage at 0.9:", free_membership_marginal["empirical_coverage"]) print("free full_cov empirical coverage at 0.9:", free_membership_full_cov["empirical_coverage"]) .. rst-class:: sphx-glr-script-out .. code-block:: none free marginal interval_type: marginal free full_cov interval_type: full_cov free marginal empirical coverage at 0.9: 0.9791666666666666 free full_cov empirical coverage at 0.9: 1.0 .. GENERATED FROM PYTHON SOURCE LINES 301-307 Plot free-orientation uncertainty representations ------------------------------------------------- The left plot shows time-aggregated component norms for one active source. The right plot compares aggregated pre-calibration curves for ``marginal`` and ``full_cov`` uncertainty. .. GENERATED FROM PYTHON SOURCE LINES 307-363 .. code-block:: Python free_source_idx = int(np.atleast_1d(active_free)[0]) true_free_norm = np.linalg.norm(np.mean(x_true_free, axis=2), axis=1) est_free_norm = np.linalg.norm(np.mean(free_result["posterior_mean_reshaped"], axis=2), axis=1) free_component_var_agg = free_membership_marginal["posterior_var_agg"][free_source_idx] fig, axes = plt.subplots(1, 2, figsize=(10, 4.2)) axes[0].bar( np.arange(3) - 0.2, np.mean(x_true_free[free_source_idx], axis=1), width=0.4, label="true aggregated", ) axes[0].bar( np.arange(3) + 0.2, np.mean(free_result["posterior_mean_reshaped"][free_source_idx], axis=1), width=0.4, label="predicted aggregated", ) axes[0].errorbar( np.arange(3) + 0.2, np.mean(free_result["posterior_mean_reshaped"][free_source_idx], axis=1), yerr=np.sqrt(free_component_var_agg), fmt="none", ecolor="black", capsize=4, label="marginal std", ) axes[0].set( xlabel="Orientation component", ylabel="Aggregated amplitude (nAm)", title=f"Free EEG: source {free_source_idx}", ) axes[0].legend(loc="best") axes[1].plot([0, 1], [0, 1], "--", color="0.5", label="perfect calibration") axes[1].plot( free_curve_marginal["nominal_coverages"], free_curve_marginal["empirical_coverages"], marker="o", label="free marginal", ) axes[1].plot( free_curve_full_cov["nominal_coverages"], free_curve_full_cov["empirical_coverages"], marker="s", label="free full_cov", ) axes[1].set( xlabel="Nominal coverage", ylabel="Empirical coverage", title="Free EEG pre-calibration curves", ) axes[1].legend(loc="best") fig.tight_layout() .. image-sg:: /auto_tutorials/images/sphx_glr_06_uncertainty_estimation_002.png :alt: Free EEG: source 12, Free EEG pre-calibration curves :srcset: /auto_tutorials/images/sphx_glr_06_uncertainty_estimation_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 364-377 Summary ------- ``UncertaintyEstimator`` is the bridge between posterior covariance and calibration-ready uncertainty. In the current workflow: - fixed orientation stores a reduced ``posterior_var`` representation; - free EEG with ``marginal`` uses pooled component-wise variances; - free EEG with ``full_cov`` uses local ``3 x 3`` covariance blocks; - aggregated calibration is the default mode, so uncertainty is evaluated on time-averaged predictions. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.536 seconds) .. _sphx_glr_download_auto_tutorials_06_uncertainty_estimation.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 06_uncertainty_estimation.ipynb <06_uncertainty_estimation.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 06_uncertainty_estimation.py <06_uncertainty_estimation.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 06_uncertainty_estimation.zip <06_uncertainty_estimation.zip>`