diff --git a/doc/tutorials/error_analysis/error_analysis_part1.ipynb b/doc/tutorials/error_analysis/error_analysis_part1.ipynb
index 8e85c210993..9251fcee29b 100644
--- a/doc/tutorials/error_analysis/error_analysis_part1.ipynb
+++ b/doc/tutorials/error_analysis/error_analysis_part1.ipynb
@@ -39,6 +39,9 @@
     "import numpy as np\n",
     "import matplotlib.pyplot as plt\n",
     "plt.rcParams.update({'font.size': 18})\n",
+    "import sys\n",
+    "import logging\n",
+    "logging.basicConfig(level=logging.INFO, stream=sys.stdout)\n",
     "\n",
     "np.random.seed(43)\n",
     "\n",
@@ -54,7 +57,7 @@
     "\n",
     "# generate simulation data using the AR(1) process\n",
     "\n",
-    "print(\"Generating data sets for the tutorial ...\\n\")\n",
+    "logging.info(\"Generating data sets for the tutorial ...\")\n",
     "\n",
     "N_SAMPLES = 100000\n",
     "\n",
@@ -69,7 +72,7 @@
     "time_series_2 = ar_1_process(N_SAMPLES, C_2, PHI_2, EPS_2)\n",
     "\n",
     "\n",
-    "print(\"Done\")"
+    "logging.info(\"Done\")"
    ]
   },
   {
@@ -292,8 +295,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "print(\"Best guess for measured quantity:\",avg)\n",
-    "print(\"Standard error of the mean:\",sem)"
+    "logging.info(f\"Best guess for measured quantity: {avg:.3f}\")\n",
+    "logging.info(f\"Standard error of the mean: {sem:.3f}\")"
    ]
   },
   {
@@ -413,7 +416,7 @@
     "plt.ylabel(\"SEM\")\n",
     "plt.show()\n",
     "\n",
-    "print(\"Final Standard Error of the Mean: {:.3f}\".format(fit_params[2]))"
+    "logging.info(f\"Final Standard Error of the Mean: {fit_params[2]:.3f}\")"
    ]
   },
   {
diff --git a/doc/tutorials/error_analysis/error_analysis_part2.ipynb b/doc/tutorials/error_analysis/error_analysis_part2.ipynb
index 4a7b097c3ea..23e789eeb86 100644
--- a/doc/tutorials/error_analysis/error_analysis_part2.ipynb
+++ b/doc/tutorials/error_analysis/error_analysis_part2.ipynb
@@ -37,6 +37,9 @@
     "import numpy as np\n",
     "import matplotlib.pyplot as plt\n",
     "plt.rcParams.update({'font.size': 18})\n",
+    "import sys\n",
+    "import logging\n",
+    "logging.basicConfig(level=logging.INFO, stream=sys.stdout)\n",
     "\n",
     "np.random.seed(43)\n",
     "\n",
@@ -52,7 +55,7 @@
     "\n",
     "# generate simulation data using the AR(1) process\n",
     "\n",
-    "print(\"Generating data sets for the tutorial ...\\n\")\n",
+    "logging.info(\"Generating data sets for the tutorial ...\\n\")\n",
     "\n",
     "N_SAMPLES = 100000\n",
     "\n",
@@ -67,7 +70,7 @@
     "time_series_2 = ar_1_process(N_SAMPLES, C_2, PHI_2, EPS_2)\n",
     "\n",
     "\n",
-    "print(\"Done\")"
+    "logging.info(\"Done\")"
    ]
   },
   {
@@ -267,7 +270,7 @@
     "plt.ylabel(\"$\\hat{R}^{XX}_j$\")\n",
     "plt.show()\n",
     "\n",
-    "print(\"Exponential autocorrelation time: {:.2f} sampling intervals\".format(popt[1]))"
+    "logging.info(f\"Exponential autocorrelation time: {popt[1]:.2f} sampling intervals\")"
    ]
   },
   {
@@ -344,7 +347,7 @@
     "plt.ylabel(r\"$\\tau_{X, \\mathrm{int}}$\")\n",
     "plt.show()\n",
     "\n",
-    "print(\"j_max = {}\".format(j_max))"
+    "logging.info(\"j_max = {}\".format(j_max))"
    ]
   },
   {
@@ -361,15 +364,15 @@
    "outputs": [],
    "source": [
     "tau_int = tau_int_v[j_max]\n",
-    "print(\"Integrated autocorrelation time: {:.2f} time steps\\n\".format(tau_int))\n",
+    "logging.info(\"Integrated autocorrelation time: {:.2f} time steps\\n\".format(tau_int))\n",
     "\n",
     "N_eff = N_SAMPLES / (2 * tau_int)\n",
-    "print(\"Effective number of samples: {:.2f}\".format(N_eff))\n",
-    "print(\"Original number of samples: {}\".format(N_SAMPLES))\n",
-    "print(\"Ratio: {:.4f}\\n\".format(N_eff/N_SAMPLES))\n",
+    "logging.info(\"Effective number of samples: {:.2f}\".format(N_eff))\n",
+    "logging.info(\"Original number of samples: {}\".format(N_SAMPLES))\n",
+    "logging.info(\"Ratio: {:.4f}\\n\".format(N_eff/N_SAMPLES))\n",
     "\n",
     "sem = np.sqrt(autocov[0]/N_eff)\n",
-    "print(\"Standard error of the mean: {:.4f}\".format(sem))"
+    "logging.info(\"Standard error of the mean: {:.4f}\".format(sem))"
    ]
   },
   {
@@ -385,7 +388,7 @@
     "  * `C` (which is the criterion to find $j_\\mathrm{max}$) and \n",
     "  * `window` (an integer that defines how much of the auto-covariance function is computed during the analysis).\n",
     "  \n",
-    "  The function shall return the SEM and print out:\n",
+    "  The function shall return the SEM and logging.info out:\n",
     "  * mean\n",
     "  * SEM\n",
     "  * integrated autocorrelation time\n",
@@ -456,11 +459,11 @@
     "    plt.title(\"\")\n",
     "    plt.show()\n",
     "    \n",
-    "    # print out stuff\n",
-    "    print(\"Mean value: {:.4f}\".format(avg))\n",
-    "    print(\"Standard error of the mean: {:.4f}\".format(sem))\n",
-    "    print(\"Integrated autocorrelation time: {:.2f} time steps\".format(tau_int))\n",
-    "    print(\"Effective number of samples: {:.2f}\".format(N_eff))\n",
+    "    # logging.info out stuff\n",
+    "    logging.info(\"Mean value: {:.4f}\".format(avg))\n",
+    "    logging.info(\"Standard error of the mean: {:.4f}\".format(sem))\n",
+    "    logging.info(\"Integrated autocorrelation time: {:.2f} time steps\".format(tau_int))\n",
+    "    logging.info(\"Effective number of samples: {:.2f}\".format(N_eff))\n",
     "\n",
     "    return sem\n",
     "\n",