diff --git a/docs/examples/emulator/robustgasp/plot_scalargasp.py b/docs/examples/emulator/robustgasp/plot_scalargasp.py
index 710b163905f8e3558e262c2969267a96c9c84051..cbb5a6598d4a31bbf5601219b92dd3ef44674e9a 100644
--- a/docs/examples/emulator/robustgasp/plot_scalargasp.py
+++ b/docs/examples/emulator/robustgasp/plot_scalargasp.py
@@ -6,7 +6,7 @@ ScalarGaSP: GP emulation for a single-output function
 
 # %% md
 # 
-# This example shows how apply Gaussian process emulation to a single-output
+# This example shows how to apply Gaussian process emulation to a single-output
 # function using class :class:`.ScalarGaSP`.
 #
 # The task is to build a GP emulator for the function :math:`y = x * sin(x)`
@@ -54,7 +54,7 @@ validation = emulator.loo_validate()
 # %% md
 #
 # Let's plot emulator predictions vs actual outputs. The error bar indicates the
-# 95% credible interval.
+# standard deviation.
 
 import matplotlib.pyplot as plt
 
@@ -66,7 +66,8 @@ ax.set_xlim(np.min(y)-1,np.max(y)+1)
 ax.set_ylim(np.min(y)-1,np.max(y)+1)
 
 _ = ax.plot([np.min(y)-1,np.max(y)+1], [np.min(y)-1,np.max(y)+1])
-_ = ax.errorbar(y, validation[:,0], validation[:,1], fmt='.', linestyle='')
+_ = ax.errorbar(y, validation[:,0], validation[:,1], fmt='.', linestyle='', label='prediction and std')
+_ = plt.legend()
 
 
 # %% md
@@ -81,7 +82,8 @@ testing_input = np.arange(0,10,0.1)
 predictions = emulator.predict(testing_input)
 
 plt.plot(testing_input, predictions[:, 0], 'r-', label= "mean")
-plt.scatter(x, y, s=15, c='k', label="training data", zorder=2)
+plt.scatter(x, y, s=15, c='k', label="training data", zorder=3)
+plt.plot(testing_input, f(testing_input), 'k:', zorder=2, label="true function")
 plt.fill_between(testing_input, predictions[:, 1], predictions[:, 2], alpha=0.3, label="95% CI")
 plt.xlabel('x')
 plt.ylabel('emulator-predicted y')
@@ -100,6 +102,7 @@ for i in range(nsamples):
     plt.plot(testing_input, samples[:,i], '--', label=f'sample{i+1}')
 
 plt.scatter(x, y, s=15, c='k', label="training data", zorder=2)
+plt.plot(testing_input, f(testing_input), 'k:', zorder=2, label="true function")
 plt.fill_between(testing_input, predictions[:, 1], predictions[:, 2], alpha=0.3, label="95% CI")
 plt.xlabel('x')
 plt.ylabel('emulator-predicted y')