diff --git a/package/CHANGELOG b/package/CHANGELOG
index 7b75455f19..1535ae56dd 100644
--- a/package/CHANGELOG
+++ b/package/CHANGELOG
@@ -21,6 +21,8 @@ The rules for this file:
  * 2.11.0
 
 Fixes
+ * Fix doctests in pca.py for newer NumPy output format and testsetup blocks
+    (Issue #3925, PR #5301)
  * Fixes TypeError with np.int64 indexing in GSD Reader (Issue #5224)
  * Fixed bug in add_transformations allowing non-callable transformations
    (Issue #2558, PR #2558)
diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py
index fddbf7d009..13d0f71a49 100644
--- a/package/MDAnalysis/analysis/pca.py
+++ b/package/MDAnalysis/analysis/pca.py
@@ -688,10 +688,10 @@ def rmsip(self, other, n_components=None):
 
         .. testsetup::
 
-            >>> import MDAnalysis as mda
-            >>> import MDAnalysis.analysis.pca as pca
-            >>> from MDAnalysis.tests.datafiles import PSF, DCD
-            >>> u = mda.Universe(PSF, DCD)
+             import MDAnalysis as mda
+             import MDAnalysis.analysis.pca as pca
+             from MDAnalysis.tests.datafiles import PSF, DCD
+             u = mda.Universe(PSF, DCD)
 
 
         You can compare the RMSIP between different intervals of the same trajectory.
@@ -703,9 +703,9 @@ def rmsip(self, other, n_components=None):
             >>> second_interval = pca.PCA(u, select="backbone").run(start=25, stop=50)
             >>> last_interval = pca.PCA(u, select="backbone").run(start=75)
             >>> round(first_interval.rmsip(second_interval, n_components=3), 6)
-            0.381476
+            np.float64(0.381476)
             >>> round(first_interval.rmsip(last_interval, n_components=3), 6)
-            0.174782
+            np.float64(0.174782)
 
 
         See also
@@ -870,10 +870,10 @@ def rmsip(a, b, n_components=None):
 
     .. testsetup::
 
-        >>> import MDAnalysis as mda
-        >>> import MDAnalysis.analysis.pca as pca
-        >>> from MDAnalysis.tests.datafiles import PSF, DCD
-        >>> u = mda.Universe(PSF, DCD)
+         import MDAnalysis as mda
+         import MDAnalysis.analysis.pca as pca
+         from MDAnalysis.tests.datafiles import PSF, DCD
+         u = mda.Universe(PSF, DCD)
 
 
     You can compare the RMSIP between different intervals of the same trajectory.
@@ -887,11 +887,11 @@ def rmsip(a, b, n_components=None):
         >>> round(pca.rmsip(first_interval.results.p_components.T,
         ...           second_interval.results.p_components.T,
         ...           n_components=3), 6)
-        0.381476
+        np.float64(0.381476)
         >>> round(pca.rmsip(first_interval.results.p_components.T,
         ...           last_interval.results.p_components.T,
         ...           n_components=3), 6)
-        0.174782
+        np.float64(0.174782)
 
 
     .. versionadded:: 1.0.0