diff --git a/demos/Wiener_Filter.ipynb b/demos/Wiener_Filter.ipynb
index c2c414ae9a81fb8935951c80daa2ab1495a95ee8..12e26c5b0e1573e11ad03a74994ea22f5d7535f9 100644
--- a/demos/Wiener_Filter.ipynb
+++ b/demos/Wiener_Filter.ipynb
@@ -93,7 +93,6 @@
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
-    "collapsed": true,
     "slideshow": {
      "slide_type": "-"
     }
@@ -161,7 +160,6 @@
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
-    "collapsed": true,
     "slideshow": {
      "slide_type": "-"
     }
@@ -225,9 +223,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {
-    "collapsed": true
-   },
+   "metadata": {},
    "outputs": [],
    "source": [
     "s_space = ift.RGSpace(N_pixels)\n",
@@ -268,7 +264,6 @@
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
-    "collapsed": true,
     "slideshow": {
      "slide_type": "-"
     }
@@ -293,7 +288,6 @@
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
-    "collapsed": true,
     "slideshow": {
      "slide_type": "-"
     }
@@ -414,7 +408,6 @@
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
-    "collapsed": true,
     "slideshow": {
      "slide_type": "skip"
     }
@@ -448,7 +441,6 @@
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
-    "collapsed": true,
     "slideshow": {
      "slide_type": "-"
     }
@@ -471,7 +463,6 @@
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
-    "collapsed": true,
     "slideshow": {
      "slide_type": "skip"
     }
@@ -503,12 +494,7 @@
    },
    "outputs": [],
    "source": [
-    "sc = ift.probing.utils.StatCalculator()\n",
-    "for i in range(200):\n",
-    "    print(i)\n",
-    "    sc.add(HT(curv.generate_posterior_sample()))\n",
-    "\n",
-    "m_var = sc.var"
+    "m_mean, m_var = ift.probe_with_posterior_samples(curv, m, HT, 200)\n"
    ]
   },
   {
@@ -526,7 +512,6 @@
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
-    "collapsed": true,
     "slideshow": {
      "slide_type": "skip"
     }
@@ -601,9 +586,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {
-    "collapsed": true
-   },
+   "metadata": {},
    "outputs": [],
    "source": [
     "N_pixels = 256      # Number of pixels\n",
@@ -661,18 +644,7 @@
     "m = D(j)\n",
     "\n",
     "# Uncertainty\n",
-    "sc = ift.probing.utils.StatCalculator()\n",
-    "\n",
-    "IC = ift.GradientNormController(iteration_limit=50000,\n",
-    "                                tol_abs_gradnorm=0.1)\n",
-    "inverter = ift.ConjugateGradient(controller=IC)\n",
-    "curv = ift.library.wiener_filter_curvature.WienerFilterCurvature(R,N,Sh,inverter)\n",
-    "\n",
-    "for i in range(20):\n",
-    "    print(i)\n",
-    "    sc.add(HT(curv.generate_posterior_sample()))\n",
-    "\n",
-    "m_var = sc.var\n",
+    "m_mean, m_var = ift.probe_with_posterior_samples(curv, m, HT, 20)\n",
     "\n",
     "# Get data\n",
     "s_power = ift.power_analyze(sh)\n",
@@ -793,21 +765,21 @@
  "metadata": {
   "celltoolbar": "Slideshow",
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Python 2",
    "language": "python",
-   "name": "python3"
+   "name": "python2"
   },
   "language_info": {
    "codemirror_mode": {
     "name": "ipython",
-    "version": 3
+    "version": 2
    },
    "file_extension": ".py",
    "mimetype": "text/x-python",
    "name": "python",
    "nbconvert_exporter": "python",
-   "pygments_lexer": "ipython3",
-   "version": "3.6.4"
+   "pygments_lexer": "ipython2",
+   "version": "2.7.12"
   }
  },
  "nbformat": 4,
diff --git a/demos/paper_demos/cartesian_wiener_filter.py b/demos/paper_demos/cartesian_wiener_filter.py
index 46e5a72881813e66e9419ca5638969bf389eb1a6..5c690f418a78b2e37aa3e62de486e537c842b661 100644
--- a/demos/paper_demos/cartesian_wiener_filter.py
+++ b/demos/paper_demos/cartesian_wiener_filter.py
@@ -108,5 +108,5 @@ if __name__ == "__main__":
     ift.plot(ift.Field(plot_space,val=data.val), name='data.png', **plotdict)
     ift.plot(ift.Field(plot_space,val=m.val), name='map.png', **plotdict)
     # sampling the uncertainty map
-    mean, variance = ift.probe_with_posterior_samples(wiener_curvature, ht, 10)
+    mean, variance = ift.probe_with_posterior_samples(wiener_curvature, m_k, ht, 10)
     ift.plot(ift.Field(plot_space, val=ift.sqrt(variance).val), name="uncertainty.png", **plotdict)
diff --git a/demos/paper_demos/wiener_filter.py b/demos/paper_demos/wiener_filter.py
index aa873d35a4fe8e311910aa002caf2e51c128af74..2d48e6c9b11d99b83436c19e62b6c9cec4f47f4c 100644
--- a/demos/paper_demos/wiener_filter.py
+++ b/demos/paper_demos/wiener_filter.py
@@ -72,7 +72,7 @@ if __name__ == "__main__":
     sample_mean = ift.Field.zeros(signal_space)
     n_samples = 10
     for i in range(n_samples):
-        sample = ht(wiener_curvature.generate_posterior_sample()) + m
+        sample = ht(wiener_curvature.draw_sample()) + m
         sample_variance += sample**2
         sample_mean += sample
     variance = sample_variance/n_samples - (sample_mean/n_samples)**2
diff --git a/nifty4/library/critical_power_energy.py b/nifty4/library/critical_power_energy.py
index 2f7f289ab0ff5b3c2891e765ce100b14d7e50ef1..bc08071314556928f6baa177aa419e3d9f46473b 100644
--- a/nifty4/library/critical_power_energy.py
+++ b/nifty4/library/critical_power_energy.py
@@ -91,7 +91,7 @@ class CriticalPowerEnergy(Energy):
             if self.D is not None:
                 w = Field.zeros(self.position.domain, dtype=self.m.dtype)
                 for i in range(self.samples):
-                    sample = self.D.generate_posterior_sample() + self.m
+                    sample = self.D.draw_sample() + self.m
                     w += P(abs(sample)**2)
 
                 w *= 1./self.samples
diff --git a/nifty4/library/noise_energy.py b/nifty4/library/noise_energy.py
index 2ea0bcb8dd3de2d296395de842f46082c09e5eaa..465b3a621d56cae0085ae1d028dff33dd618b874 100644
--- a/nifty4/library/noise_energy.py
+++ b/nifty4/library/noise_energy.py
@@ -46,7 +46,7 @@ class NoiseEnergy(Energy):
             if samples is None or samples == 0:
                 xi_sample_list = [xi]
             else:
-                xi_sample_list = [D.generate_posterior_sample() + xi
+                xi_sample_list = [D.draw_sample() + xi
                                   for _ in range(samples)]
         self.xi_sample_list = xi_sample_list
         self.inverter = inverter
diff --git a/nifty4/library/nonlinear_power_energy.py b/nifty4/library/nonlinear_power_energy.py
index c7a0a4616cda1f54108997d656d9ef7747fe9a0f..670757d1f64c25efa18dc471e02eafe41c3b2034 100644
--- a/nifty4/library/nonlinear_power_energy.py
+++ b/nifty4/library/nonlinear_power_energy.py
@@ -70,7 +70,7 @@ class NonlinearPowerEnergy(Energy):
             if samples is None or samples == 0:
                 xi_sample_list = [xi]
             else:
-                xi_sample_list = [D.generate_posterior_sample() + xi
+                xi_sample_list = [D.draw_sample() + xi
                                   for _ in range(samples)]
         self.xi_sample_list = xi_sample_list
         self.inverter = inverter
diff --git a/nifty4/library/wiener_filter_curvature.py b/nifty4/library/wiener_filter_curvature.py
index 53f0c708250944ec1c659149558b455d21b7ccd2..25f1f342f13ce661681434b46e1582c7d42bb603 100644
--- a/nifty4/library/wiener_filter_curvature.py
+++ b/nifty4/library/wiener_filter_curvature.py
@@ -59,23 +59,22 @@ class WienerFilterCurvature(EndomorphicOperator):
     def apply(self, x, mode):
         return self._op.apply(x, mode)
 
-    def generate_posterior_sample(self):
-        """ Generates a posterior sample from a Gaussian distribution with
-        given mean and covariance.
+    def draw_sample(self):
+        """ Generates a sample from a Gaussian distribution with
+        covariance given by the operator.
 
         This method generates samples by setting up the observation and
         reconstruction of a mock signal in order to obtain residuals of the
-        right correlation which are added to the given mean.
+        right correlation.
 
         Returns
         -------
         sample : Field
-            Returns the a sample from the Gaussian of mean zero and
-            given covariance.
+            Returns the a sample from the Gaussian of given covariance.
         """
 
-        mock_signal = self.S.generate_posterior_sample()
-        mock_noise = self.N.generate_posterior_sample()
+        mock_signal = self.S.draw_sample()
+        mock_noise = self.N.draw_sample()
 
         mock_data = self.R(mock_signal) + mock_noise
 
diff --git a/nifty4/operators/diagonal_operator.py b/nifty4/operators/diagonal_operator.py
index 8be3c7a6805c746cfb3f0650b13db74c1a500bc3..15c1418d432c8dc6b86a9fe3dab13d55741d3185 100644
--- a/nifty4/operators/diagonal_operator.py
+++ b/nifty4/operators/diagonal_operator.py
@@ -141,19 +141,18 @@ class DiagonalOperator(EndomorphicOperator):
         return DiagonalOperator(self._diagonal.conjugate(), self._domain,
                                 self._spaces)
 
-    def generate_posterior_sample(self):
-        """ Generates a posterior sample from a Gaussian distribution with
-        given mean and covariance.
+    def draw_sample(self):
+        """ Generates a sample from a Gaussian distribution with
+        covariance given by the operator.
 
         This method generates samples by setting up the observation and
         reconstruction of a mock signal in order to obtain residuals of the
-        right correlation which are added to the given mean.
+        right correlation.
 
         Returns
         -------
         sample : Field
-            Returns the a sample from the Gaussian of given mean and
-            covariance.
+            Returns the a sample from the Gaussian of given covariance.
         """
 
         if self._spaces is not None:
diff --git a/nifty4/operators/scaling_operator.py b/nifty4/operators/scaling_operator.py
index 678788e75614ecaa987ce0459111b9ce326867a9..e3429048e7f28b9a4537bb91092c476bc65cbdd5 100644
--- a/nifty4/operators/scaling_operator.py
+++ b/nifty4/operators/scaling_operator.py
@@ -89,19 +89,18 @@ class ScalingOperator(EndomorphicOperator):
         return (self.TIMES | self.ADJOINT_TIMES |
                 self.INVERSE_TIMES | self.ADJOINT_INVERSE_TIMES)
 
-    def generate_posterior_sample(self):
-        """ Generates a posterior sample from a Gaussian distribution with
-        given mean and covariance.
+    def draw_sample(self):
+        """ Generates a sample from a Gaussian distribution with
+        covariance given by the operator.
 
         This method generates samples by setting up the observation and
         reconstruction of a mock signal in order to obtain residuals of the
-        right correlation which are added to the given mean.
+        right correlation.
 
         Returns
         -------
         sample : Field
-            Returns the a sample from the Gaussian of given mean and
-            covariance.
+            Returns the a sample from the Gaussian of given covariance.
         """
 
         return Field.from_random(random_type="normal",
diff --git a/nifty4/probing/utils.py b/nifty4/probing/utils.py
index a179e55dc1ad9998558e2ca18c1ef861636656af..7aa959666128273016b642cb22649ed0ea008ff1 100644
--- a/nifty4/probing/utils.py
+++ b/nifty4/probing/utils.py
@@ -50,7 +50,7 @@ class StatCalculator(object):
 def probe_with_posterior_samples(op, m, post_op, nprobes):
     sc = StatCalculator()
     for i in range(nprobes):
-        sample = post_op(op.generate_posterior_sample() + m)
+        sample = post_op(op.draw_sample() + m)
         sc.add(sample)
 
     if nprobes == 1: