diff --git a/nimare/meta/cbmr.py b/nimare/meta/cbmr.py
index 370d31bde..156914f81 100644
--- a/nimare/meta/cbmr.py
+++ b/nimare/meta/cbmr.py
@@ -20,7 +20,13 @@
 from nimare.diagnostics import FocusFilter
 from nimare.estimator import Estimator
 from nimare.meta import models
-from nimare.utils import b_spline_bases, dummy_encoding_moderators, get_masker, mm2vox
+from nimare.utils import (
+    b_spline_bases,
+    dummy_encoding_moderators,
+    get_masker,
+    mm2vox,
+    robust_inverse,
+)
 
 LGR = logging.getLogger(__name__)
 __version__ = _version.get_versions()["version"]
@@ -106,6 +112,8 @@ class CBMREstimator(Estimator):
 
     Notes
     -----
+    The CBMR framework was developed in :footcite:t:`yu2024neuroimaging`,
+
     Available correction methods: :meth:`~nimare.meta.cbmr.CBMRInference`.
     """
 
@@ -751,7 +759,7 @@ def _glh_con_group(self):
                 self.estimator.inputs_["foci_per_voxel"],
                 self.estimator.inputs_["foci_per_study"],
             )
-            cov_spatial_coef = np.linalg.inv(f_spatial_coef)
+            cov_spatial_coef = robust_inverse(f_spatial_coef)
             # compute numerator: contrast vector * group-wise log spatial intensity
             involved_log_intensity_per_voxel = list()
             for group in con_group_involved:
@@ -925,8 +933,7 @@ def _glh_con_moderator(self):
                 self.estimator.inputs_["foci_per_voxel"],
                 self.estimator.inputs_["foci_per_study"],
             )
-
-            cov_moderator_coef = np.linalg.inv(f_moderator_coef)
+            cov_moderator_coef = robust_inverse(f_moderator_coef)
             if m_con_moderator == 1:  # a single contrast vector, use Wald test
                 var_moderator_coef = np.diag(cov_moderator_coef)
                 involved_var_moderator_coef = con_moderator**2 @ var_moderator_coef
diff --git a/nimare/tests/conftest.py b/nimare/tests/conftest.py
index 93a539f72..d33e7179c 100644
--- a/nimare/tests/conftest.py
+++ b/nimare/tests/conftest.py
@@ -77,7 +77,7 @@ def testdata_cbmr_simulated():
     """Simulate coordinate-based dataset for tests."""
     # simulate
     ground_truth_foci, dset = create_coordinate_dataset(
-        foci=10, sample_size=(20, 40), n_studies=1000, seed=100
+        foci=10, sample_size=(20, 40), n_studies=500, seed=100
     )
     # set up group columns: diagnosis & drug_status
     n_rows = dset.annotations.shape[0]
diff --git a/nimare/tests/test_meta_cbmr.py b/nimare/tests/test_meta_cbmr.py
index 8fb243dbb..127355dbe 100644
--- a/nimare/tests/test_meta_cbmr.py
+++ b/nimare/tests/test_meta_cbmr.py
@@ -262,6 +262,7 @@ def test_StandardizeField(testdata_cbmr_simulated):
 
 
 @pytest.mark.cbmr_importerror
+@pytest.mark.skipif(TORCH_INSTALLED, reason="Torch is installed, ImportError test not applicable")
 def test_cbmr_importerror():
     """Test that ImportErrors are raised when torch is not installed."""
     with pytest.raises(ImportError):
diff --git a/nimare/utils.py b/nimare/utils.py
index 3878fadfb..eece72da7 100755
--- a/nimare/utils.py
+++ b/nimare/utils.py
@@ -207,7 +207,7 @@ def mm2vox(xyz, affine):
     From here:
     http://blog.chrisgorgolewski.org/2014/12/how-to-convert-between-voxel-and-mm.html
     """
-    with np.errstate(invalid='ignore'):
+    with np.errstate(invalid="ignore"):
         ijk = nib.affines.apply_affine(np.linalg.inv(affine), xyz).astype(int)
     return ijk
 
@@ -1366,3 +1366,39 @@ def dummy_encoding_moderators(dataset_annotations, moderators):
         else:
             new_moderators.append(moderator)
     return dataset_annotations, new_moderators
+
+
+def robust_inverse(FI, eps=1e-8):
+    """
+    Compute the robust inverse of a Fisher information matrix.
+
+    This function computes the inverse of a Fisher information matrix (FI) using
+    Singular Value Decomposition (SVD) and a thresholding approach to handle
+    small singular values, which can lead to numerical instability.
+
+    Parameters
+    ----------
+    FI : :obj:`numpy.ndarray`
+        The Fisher information matrix to be inverted.
+    eps : :obj:`float`, optional
+        A small value to threshold the singular values. Singular values below this
+        threshold will be set to zero in the inverse computation to avoid numerical
+        instability. Default is 1e-8.
+
+    Returns
+    -------
+    FI_inv : :obj:`numpy.ndarray`
+        The robust inverse of the Fisher information matrix.
+    """
+    # Improve numerical stability of FI by symmetrizing it
+    FI = (FI + FI.T) / 2
+    # Perform Singular Value Decomposition
+    # and compute the inverse using a threshold on singular values
+    # to avoid division by zero or very small values
+    U, S, VT = np.linalg.svd(FI, full_matrices=False)
+    M = S > eps
+    S_inv = S**-1
+    U = ((U + VT.T) / 2) * M[None, :]
+    # Set small singular values to zero in the inverse
+    FI_inv = U @ np.diag(S_inv) @ U.T
+    return FI_inv