diff --git a/demos/Wiener_Filter.ipynb b/demos/Wiener_Filter.ipynb
index 12e26c5b0e1573e11ad03a74994ea22f5d7535f9..c576c81a0745e166494f6566fa4b2054b39c8481 100644
--- a/demos/Wiener_Filter.ipynb
+++ b/demos/Wiener_Filter.ipynb
@@ -80,8 +80,8 @@
    "source": [
     "## Wiener Filter: Example\n",
     "\n",
-    "- One-dimensional signal with power spectrum: $$P(k) = P_0\\,\\left(1+\\left(\\frac{k}{k_0}\\right)^2\\right)^{-\\gamma /2},$$\n",
-    "with $P_0 = 0.2, k_0 = 5, \\gamma = 4$. Recall: $P(k)$ defines an isotropic and homogeneous $S$.\n",
+    "- We assume statistical homogeneity and isotropy. Therefore the signal covariance $S$ is diagonal in harmonic space, and is described by a one-dimensional power spectrum, assumed here as $$P(k) = P_0\\,\\left(1+\\left(\\frac{k}{k_0}\\right)^2\\right)^{-\\gamma /2},$$\n",
+    "with $P_0 = 0.2, k_0 = 5, \\gamma = 4$.\n",
     "- $N = 0.2 \\cdot \\mathbb{1}$.\n",
     "- Number of data points $N_{pix} = 512$.\n",
     "- reconstruction in harmonic space.\n",
@@ -328,9 +328,9 @@
    "outputs": [],
    "source": [
     "plt.figure(figsize=(15,10))\n",
-    "plt.plot(s_data, 'g', label=\"Signal\")\n",
-    "plt.plot(d_data, 'k+', label=\"Data\")\n",
-    "plt.plot(m_data, 'r', label=\"Reconstruction\")\n",
+    "plt.plot(s_data, 'r', label=\"Signal\", linewidth=3)\n",
+    "plt.plot(d_data, 'k.', label=\"Data\")\n",
+    "plt.plot(m_data, 'k', label=\"Reconstruction\",linewidth=3)\n",
     "plt.title(\"Reconstruction\")\n",
     "plt.legend()\n",
     "plt.show()"
@@ -347,9 +347,9 @@
    "outputs": [],
    "source": [
     "plt.figure(figsize=(15,10))\n",
-    "plt.plot(s_data - s_data, 'g', label=\"Signal\")\n",
-    "plt.plot(d_data - s_data, 'k+', label=\"Data\")\n",
-    "plt.plot(m_data - s_data, 'r', label=\"Reconstruction\")\n",
+    "plt.plot(s_data - s_data, 'r', label=\"Signal\", linewidth=3)\n",
+    "plt.plot(d_data - s_data, 'k.', label=\"Data\")\n",
+    "plt.plot(m_data - s_data, 'k', label=\"Reconstruction\",linewidth=3)\n",
     "plt.axhspan(-noise_amplitude,noise_amplitude, facecolor='0.9', alpha=.5)\n",
     "plt.title(\"Residuals\")\n",
     "plt.legend()\n",
@@ -383,9 +383,9 @@
     "ymin = min(m_power_data)\n",
     "plt.ylim(ymin, 1)\n",
     "xs = np.arange(1,int(N_pixels/2),.1)\n",
-    "plt.plot(xs, pow_spec(xs), label=\"True Power Spectrum\", linewidth=.7, color='k')\n",
-    "plt.plot(s_power_data, 'g', label=\"Signal\")\n",
-    "plt.plot(m_power_data, 'r', label=\"Reconstruction\")\n",
+    "plt.plot(xs, pow_spec(xs), label=\"True Power Spectrum\", color='k',alpha=0.5)\n",
+    "plt.plot(s_power_data, 'r', label=\"Signal\")\n",
+    "plt.plot(m_power_data, 'k', label=\"Reconstruction\")\n",
     "plt.axhline(noise_amplitude**2 / N_pixels, color=\"k\", linestyle='--', label=\"Noise level\", alpha=.5)\n",
     "plt.axhspan(noise_amplitude**2 / N_pixels, ymin, facecolor='0.9', alpha=.5)\n",
     "plt.title(\"Power Spectrum\")\n",
@@ -545,8 +545,8 @@
    "outputs": [],
    "source": [
     "plt.figure(figsize=(15,10))\n",
-    "plt.plot(s_data, 'g', label=\"Signal\", linewidth=1)\n",
-    "plt.plot(d_data, 'k+', label=\"Data\", alpha=1)\n",
+    "plt.plot(s_data, 'r', label=\"Signal\", linewidth=3)\n",
+    "plt.plot(d_data, 'k.', label=\"Data\")\n",
     "plt.axvspan(l, h, facecolor='0.8', alpha=.5)\n",
     "plt.title(\"Incomplete Data\")\n",
     "plt.legend()"
@@ -563,11 +563,11 @@
    "outputs": [],
    "source": [
     "fig = plt.figure(figsize=(15,10))\n",
-    "plt.plot(s_data, 'g', label=\"Signal\", alpha=1, linewidth=4)\n",
-    "plt.plot(d_data, 'k+', label=\"Data\", alpha=.5)\n",
-    "plt.plot(m_data, 'r', label=\"Reconstruction\")\n",
-    "plt.axvspan(l, h, facecolor='0.8', alpha=.5)\n",
-    "plt.fill_between(range(N_pixels), m_data - uncertainty, m_data + uncertainty, facecolor='0')\n",
+    "plt.axvspan(l, h, facecolor='0.8',alpha=0.5)\n",
+    "plt.fill_between(range(N_pixels), m_data - uncertainty, m_data + uncertainty, facecolor='0.5', alpha=0.5)\n",
+    "plt.plot(s_data, 'r', label=\"Signal\", alpha=1, linewidth=3)\n",
+    "plt.plot(d_data, 'k.', label=\"Data\")\n",
+    "plt.plot(m_data, 'k', label=\"Reconstruction\", linewidth=3)\n",
     "plt.title(\"Reconstruction of incomplete data\")\n",
     "plt.legend()"
    ]
@@ -702,10 +702,12 @@
     "mi = np.min(s_data)\n",
     "ma = np.max(s_data)\n",
     "\n",
-    "fig, axes = plt.subplots(2, 2, figsize=(15, 15))\n",
+    "fig, axes = plt.subplots(3, 2, figsize=(15, 22.5))\n",
+    "samp1 = HT(curv.draw_sample()+m).val\n",
+    "samp2 = HT(curv.draw_sample()+m).val\n",
     "\n",
-    "data = [s_data, m_data, s_data - m_data, uncertainty]\n",
-    "caption = [\"Signal\", \"Reconstruction\", \"Residuals\", \"Uncertainty Map\"]\n",
+    "data = [s_data, m_data, samp1, samp2, s_data - m_data, uncertainty]\n",
+    "caption = [\"Signal\", \"Reconstruction\", \"Sample 1\", \"Sample 2\", \"Residuals\", \"Uncertainty Map\"]\n",
     "\n",
     "for ax in axes.flat:\n",
     "    im = ax.imshow(data.pop(0), interpolation='nearest', cmap=cm, vmin=mi, vmax=ma)\n",