diff --git a/estimation/estimation_dynamical_models.ipynb b/estimation/estimation_dynamical_models.ipynb
index 1d79942..e1911a4 100644
--- a/estimation/estimation_dynamical_models.ipynb
+++ b/estimation/estimation_dynamical_models.ipynb
@@ -647,7 +647,7 @@
     "ax1.plot(time2plt_normalized, (mex_prop[:, 2] - mex_sim_obs[:, 2]), label=r'$\\Delta z$')\n",
     "ax1.plot(time2plt_normalized, np.linalg.norm((mex_prop[:, 0:3] - mex_sim_obs[:, 0:3]), axis=1), label=r'$||\\Delta X||$')\n",
     "\n",
-    "ax1.set_title(\"Element-wise difference with propagated observations\")\n",
+    "ax1.set_title(\"Element-wise difference between true and estimated states\")\n",
     "ax1.set_xlabel(r'$Time$ [days]')\n",
     "ax1.set_ylabel(r'$\\Delta X$ [m]')\n",
     "ax1.legend()\n",
diff --git a/estimation/estimation_dynamical_models.py b/estimation/estimation_dynamical_models.py
index eeef12d..68b9c5f 100644
--- a/estimation/estimation_dynamical_models.py
+++ b/estimation/estimation_dynamical_models.py
@@ -122,9 +122,9 @@
 
 ## Set Up the Observations
 """
+Having set the underlying environment model of the simulated orbit, we can define the observational model. This entails the addition all required ground stations, the definition of the observation links and types, as well as the precise simulation settings.
 """
 
-Having set the underlying environment model of the simulated orbit, we can define the observational model. This entails the addition all required ground stations, the definition of the observation links and types, as well as the precise simulation settings.
 ### Add a ground station
 """
 Following its real-world counterpart, our simulated Mars Express spacecraft will also be tracked using ESA's New Norcia (NNO) ESTRACK ground station. Located in North-East Australia, it will be set up with an altitude of 252m, 31.0482°S, 116.191°E.
@@ -433,7 +433,7 @@
 ax1.plot(time2plt_normalized, (mex_prop[:, 2] - mex_sim_obs[:, 2]), label=r'$\Delta z$')
 ax1.plot(time2plt_normalized, np.linalg.norm((mex_prop[:, 0:3] - mex_sim_obs[:, 0:3]), axis=1), label=r'$||\Delta X||$')
 
-ax1.set_title("Element-wise difference with propagated observations")
+ax1.set_title("Element-wise difference between true and estimated states")
 ax1.set_xlabel(r'$Time$ [days]')
 ax1.set_ylabel(r'$\Delta X$ [m]')
 ax1.legend()
diff --git a/estimation/full_estimation_example.ipynb b/estimation/full_estimation_example.ipynb
index 2a75321..ce89b1c 100644
--- a/estimation/full_estimation_example.ipynb
+++ b/estimation/full_estimation_example.ipynb
@@ -33,9 +33,11 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": null,
    "id": "b71677c3",
-   "metadata": {},
+   "metadata": {
+    "is_executing": true
+   },
    "outputs": [],
    "source": [
     "# Load required standard modules\n",
@@ -538,6 +540,13 @@
     "# Save the true parameters to later analyse the error\n",
     "truth_parameters = parameters_to_estimate.parameter_vector\n",
     "\n",
+    "# Perturb the initial state estimate from the truth (10 m in position; 0.1 m/s in velocity)\n",
+    "perturbed_parameters = truth_parameters.copy( )\n",
+    "for i in range(3):\n",
+    "    perturbed_parameters[i] += 10.0\n",
+    "    perturbed_parameters[i+3] += 0.01\n",
+    "parameters_to_estimate.parameter_vector = perturbed_parameters\n",
+    "\n",
     "# Create input object for the estimation\n",
     "convergence_checker = estimation.estimation_convergence_checker(maximum_iterations=4)\n",
     "estimation_input = estimation.EstimationInput(\n",
@@ -658,7 +667,7 @@
     "\n",
     "plt.figure(figsize=(9, 5))\n",
     "plt.title(\"Observations as a function of time\")\n",
-    "plt.scatter(observation_times / 3600.0, observations_list * constants.SPEED_OF_LIGHT)\n",
+    "plt.scatter(observation_times / 3600.0, observations_list)\n",
     "\n",
     "plt.xlabel(\"Time [hr]\")\n",
     "plt.ylabel(\"Range rate [m/s]\")\n",
@@ -697,16 +706,18 @@
    "source": [
     "residual_history = estimation_output.residual_history\n",
     "\n",
-    "fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(9, 6), sharex=True, sharey=True)\n",
+    "fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(9, 6))\n",
     "subplots_list = [ax1, ax2, ax3, ax4]\n",
     "\n",
     "for i in range(4):\n",
     "    subplots_list[i].scatter(observation_times, residual_history[:, i])\n",
+    "    subplots_list[i].set_ylabel(\"Observation Residual [m/s]\")\n",
+    "    subplots_list[i].set_title(\"Iteration \"+str(i+1))\n",
+    "\n",
     "\n",
     "ax3.set_xlabel(\"Time since J2000 [s]\")\n",
     "ax4.set_xlabel(\"Time since J2000 [s]\")\n",
-    "ax1.set_ylabel(\"Observation Residual [m/s]\")\n",
-    "ax3.set_ylabel(\"Observation Residual [m/s]\")\n",
+    "\n",
     "\n",
     "plt.tight_layout()\n",
     "plt.show()"
@@ -762,6 +773,9 @@
     "\n",
     "plt.figure(figsize=(9,5))\n",
     "plt.hist(final_residuals, 25)\n",
+    "plt.xlabel('Final iteration range-rate residual [m/s]')\n",
+    "plt.ylabel('Occurences [-]')\n",
+    "plt.title('Histogram of residuals on final iteration')\n",
     "\n",
     "plt.tight_layout()\n",
     "plt.show()"
diff --git a/estimation/full_estimation_example.py b/estimation/full_estimation_example.py
index a564bc7..fdb0efa 100644
--- a/estimation/full_estimation_example.py
+++ b/estimation/full_estimation_example.py
@@ -365,6 +365,14 @@
 # Save the true parameters to later analyse the error
 truth_parameters = parameters_to_estimate.parameter_vector
 
+# Perturb the initial state estimate from the truth (10 m in position; 0.1 m/s in velocity)
+perturbed_parameters = truth_parameters.copy( )
+for i in range(3):
+    perturbed_parameters[i] += 10.0
+    perturbed_parameters[i+3] += 0.01
+parameters_to_estimate.parameter_vector = perturbed_parameters
+
+
 # Create input object for the estimation
 convergence_checker = estimation.estimation_convergence_checker(maximum_iterations=4)
 estimation_input = estimation.EstimationInput(
@@ -412,7 +420,7 @@
 
 plt.figure(figsize=(9, 5))
 plt.title("Observations as a function of time")
-plt.scatter(observation_times / 3600.0, observations_list * constants.SPEED_OF_LIGHT)
+plt.scatter(observation_times / 3600.0, observations_list )
 
 plt.xlabel("Time [hr]")
 plt.ylabel("Range rate [m/s]")
@@ -429,16 +437,16 @@
 
 residual_history = estimation_output.residual_history
 
-fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(9, 6), sharex=True, sharey=True)
+fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(9, 6), sharex=True)
 subplots_list = [ax1, ax2, ax3, ax4]
 
 for i in range(4):
     subplots_list[i].scatter(observation_times, residual_history[:, i])
+    subplots_list[i].set_ylabel("Observation Residual [m/s]")
+    subplots_list[i].set_title("Iteration "+str(i+1))
 
 ax3.set_xlabel("Time since J2000 [s]")
 ax4.set_xlabel("Time since J2000 [s]")
-ax1.set_ylabel("Observation Residual [m/s]")
-ax3.set_ylabel("Observation Residual [m/s]")
 
 plt.tight_layout()
 plt.show()
@@ -459,7 +467,9 @@
 
 plt.figure(figsize=(9,5))
 plt.hist(final_residuals, 25)
-
+plt.xlabel('Final iteration range-rate residual [m/s]')
+plt.ylabel('Occurences [-]')
+plt.title('Histogram of residuals on final iteration')
 plt.tight_layout()
 plt.show()