diff --git a/demos/multi_amplitudes_consistency.py b/demos/multi_amplitudes_consistency.py
index f9ed008631f627e69a08b6a9f822acfbe1e6d586..e7b4450819a4e55099aedd6ce183e959ad470a4c 100644
--- a/demos/multi_amplitudes_consistency.py
+++ b/demos/multi_amplitudes_consistency.py
@@ -5,6 +5,11 @@ import numpy as np
 
 
 def testAmplitudesConsistency(seed, sspace):
+    def stats(op,samples):
+        sc = ift.StatCalculator()
+        for s in samples:
+            sc.add(op(s.extract(op.domain)))
+        return sc.mean.to_global_data(), sc.var.sqrt().to_global_data()
     np.random.seed(seed)
     offset_std = 30
     intergated_fluct_std0 = .003
@@ -23,13 +28,13 @@ def testAmplitudesConsistency(seed, sspace):
     op = fa.finalize(offset_std, 1E-8, '')
 
     samples = [ift.from_random('normal',op.domain) for _ in range(nsam)]
-    tot_flm, _ = fa.stats(fa.total_fluctuation,samples)
-    offset_std,_ = fa.stats(fa.amplitude_total_offset,samples)
-    intergated_fluct_std0,_ = fa.stats(fa.average_fluctuation(0),samples)
-    intergated_fluct_std1,_ = fa.stats(fa.average_fluctuation(1),samples)
+    tot_flm, _ = stats(fa.total_fluctuation,samples)
+    offset_std,_ = stats(fa.amplitude_total_offset,samples)
+    intergated_fluct_std0,_ = stats(fa.average_fluctuation(0),samples)
+    intergated_fluct_std1,_ = stats(fa.average_fluctuation(1),samples)
     
-    slice_fluct_std0,_ = fa.stats(fa.slice_fluctuation(0),samples)
-    slice_fluct_std1,_ = fa.stats(fa.slice_fluctuation(1),samples)
+    slice_fluct_std0,_ = stats(fa.slice_fluctuation(0),samples)
+    slice_fluct_std1,_ = stats(fa.slice_fluctuation(1),samples)
 
     sams = [op(s) for s in samples]
     fluct_total = fa.total_fluctuation_realized(sams)
@@ -76,7 +81,7 @@ def testAmplitudesConsistency(seed, sspace):
     fa.add_fluctuations(sspace, x, 1.5, 1.1, 2., 2.1, .5,
                         -2, 1., 'spatial', 0)
     op = fa.finalize(offset_std, .1, '')
-    em, estd = fa.stats(fa.slice_fluctuation(0),samples)
+    em, estd = stats(fa.slice_fluctuation(0),samples)
     print("Forced   slice fluct. space Std: "+str(m))
     print("Expected slice fluct. Std: " + str(em))
     np.testing.assert_allclose(m, em, rtol=0.5)
diff --git a/nifty5/library/correlated_fields.py b/nifty5/library/correlated_fields.py
index 7940747be8afc0c182c61cc6cbc34b0a6e3887d6..6659f822b862848b2b9ed65ee03f1fd9fa134a3f 100644
--- a/nifty5/library/correlated_fields.py
+++ b/nifty5/library/correlated_fields.py
@@ -444,7 +444,7 @@ class CorrelatedFieldMaker:
     def offset_amplitude_realized(samples):
         res = 0.
         for s in samples:
-            res += s.mean()**2
+            res = res + s.mean()**2
         return np.sqrt(res/len(samples))
 
     @staticmethod
@@ -461,9 +461,13 @@ class CorrelatedFieldMaker:
             return _total_fluctuation_realized(samples)
         res1, res2 = 0., 0.
         for s in samples:
-            res1 += s**2
-            res2 += s.mean(space)**2
-        return np.sqrt((res1 - res2).mean()/len(samples))
+            res1 = res1 + s**2
+            res2 = res2 + s.mean(space)**2
+        res1 = res1/len(samples)
+        res2 = res2/len(samples)
+        res = res1.mean() - res2.mean()
+        return np.sqrt(res)
+
 
     @staticmethod
     def average_fluctuation_realized(samples, space):
diff --git a/test/test_operators/test_correlated_fields.py b/test/test_operators/test_correlated_fields.py
index 6ee24f31f15a72e5d05c41d969a27a67053f555a..4beb608244b1791296611f6b8b19a4334177b8ab 100644
--- a/test/test_operators/test_correlated_fields.py
+++ b/test/test_operators/test_correlated_fields.py
@@ -32,6 +32,11 @@ import nifty5 as ift
 @pytest.mark.parametrize('Astds', [[1.,3.],[0.2,1.4]])
 @pytest.mark.parametrize('offset_std', [1.,10.])
 def testAmplitudesConsistency(rseed, sspace, Astds, offset_std):
+    def stats(op,samples):
+        sc = ift.StatCalculator()
+        for s in samples:
+            sc.add(op(s.extract(op.domain)))
+        return sc.mean.to_global_data(), sc.var.sqrt().to_global_data()
     seed(rseed)
     nsam = 100
 
@@ -45,13 +50,13 @@ def testAmplitudesConsistency(rseed, sspace, Astds, offset_std):
     op = fa.finalize(offset_std, 1E-8, '')
 
     samples = [ift.from_random('normal',op.domain) for _ in range(nsam)]
-    tot_flm, _ = fa.stats(fa.total_fluctuation,samples)
-    offset_std,_ = fa.stats(fa.amplitude_total_offset,samples)
-    intergated_fluct_std0,_ = fa.stats(fa.average_fluctuation(0),samples)
-    intergated_fluct_std1,_ = fa.stats(fa.average_fluctuation(1),samples)
+    tot_flm, _ = stats(fa.total_fluctuation,samples)
+    offset_std,_ = stats(fa.amplitude_total_offset,samples)
+    intergated_fluct_std0,_ = stats(fa.average_fluctuation(0),samples)
+    intergated_fluct_std1,_ = stats(fa.average_fluctuation(1),samples)
     
-    slice_fluct_std0,_ = fa.stats(fa.slice_fluctuation(0),samples)
-    slice_fluct_std1,_ = fa.stats(fa.slice_fluctuation(1),samples)
+    slice_fluct_std0,_ = stats(fa.slice_fluctuation(0),samples)
+    slice_fluct_std1,_ = stats(fa.slice_fluctuation(1),samples)
 
     sams = [op(s) for s in samples]
     fluct_total = fa.total_fluctuation_realized(sams)
@@ -76,7 +81,7 @@ def testAmplitudesConsistency(rseed, sspace, Astds, offset_std):
     fa.add_fluctuations(sspace, x, 1.5, 1.1, 2., 2.1, .5,
                         -2, 1., 'spatial', 0)
     op = fa.finalize(offset_std, .1, '')
-    em, estd = fa.stats(fa.slice_fluctuation(0),samples)
+    em, estd = stats(fa.slice_fluctuation(0),samples)
 
     assert_allclose(m, em, rtol=0.5)