diff --git a/docs/examples/emulator/README.rst b/docs/examples/emulator/README.rst
deleted file mode 100644
index 2665689ec0f8d8d780ca2c436c3832902722f986..0000000000000000000000000000000000000000
--- a/docs/examples/emulator/README.rst
+++ /dev/null
@@ -1,2 +0,0 @@
-emulator Examples
-^^^^^^^^^^^^^^^^^
\ No newline at end of file
diff --git a/docs/examples/emulator/robustgasp/README.rst b/docs/examples/emulator/robustgasp/README.rst
new file mode 100644
index 0000000000000000000000000000000000000000..20c048018c4ed9ea4b82b6534cabd70cb13a836c
--- /dev/null
+++ b/docs/examples/emulator/robustgasp/README.rst
@@ -0,0 +1,2 @@
+Emulator - RobustGaSP
+---------------------
\ No newline at end of file
diff --git a/docs/examples/emulator/robustgasp/plot_scalargasp.py b/docs/examples/emulator/robustgasp/plot_scalargasp.py
new file mode 100644
index 0000000000000000000000000000000000000000..e4143bc1a667a01aa47d4038aec3252183843bd8
--- /dev/null
+++ b/docs/examples/emulator/robustgasp/plot_scalargasp.py
@@ -0,0 +1,107 @@
+"""
+
+ScalarGaSP: GP emulation for a single-output function
+=====================================================
+"""
+
+# %% md
+# 
+# This example shows how 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)`
+# based on a few number of training data.
+
+# %% md
+# 
+# First, import the class :class:`.ScalarGaSP` by
+
+from psimpy.emulator import ScalarGaSP
+
+# %% md
+# 
+# Then, create an instance of :class:`.ScalarGaSP`.  The parameter ``ndim``
+# (dimension of function input ``x``) must be specified. Optional parameters, such
+# as ``method``, ``kernel_type``, etc., can be set up if desired. Here, we leave
+# all the optional parameters to their default values.
+
+emulator = ScalarGaSP(ndim=1)
+
+# %% md
+# 
+# Given training input points ``design`` and corresponding output values ``response``,
+# the emulator can be trained. Below we train an emulator using :math:`8` selected
+# points.
+
+import numpy as np
+
+def f(x):
+    #return x + 3*np.sin(x/2)
+    return x*np.sin(x)
+
+x = np.arange(2,10,1)
+y = f(x)
+
+emulator.train(design=x, response=y)
+
+# %% md
+# 
+# We can validate the performance of the trained emulator using the leave-one-out
+# cross validation method :py:meth:`loo_validate()`.
+
+validation = emulator.loo_validate()
+
+# %% md
+#
+# Let's plot emulator predictions vs actual outputs. The error bar indicates the
+# 95% credible interval.
+
+import matplotlib.pyplot as plt
+
+fig , ax = plt.subplots(figsize=(4,4))
+
+ax.set_xlabel('Actual y')
+ax.set_ylabel('Emulator-predicted y')
+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='')
+
+
+# %% md
+#
+# With the trained emulator at our deposit, we can use the
+# :py:meth:`.ScalarGaSP.predict()` method to make predictions at 
+# any arbitrary set of input points (``testing_input``). It should be noted that,
+# ``testing_trend`` should be set according to ``trend`` used during emulator
+# training. 
+
+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.fill_between(testing_input, predictions[:, 1], predictions[:, 2], alpha=0.3, label="95% CI")
+plt.xlabel('x')
+plt.ylabel('emulator-predicted y')
+plt.xlim(testing_input[0], testing_input[-1])
+_ = plt.legend()
+
+# %% md
+#
+# We can also draw any number of samples at ``testing_input``` using
+# :py:meth:`.ScalarGaSP.sample()` method.
+
+nsamples = 5
+samples = emulator.sample(testing_input, nsamples=nsamples)
+
+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.fill_between(testing_input, predictions[:, 1], predictions[:, 2], alpha=0.3, label="95% CI")
+plt.xlabel('x')
+plt.ylabel('emulator-predicted y')
+plt.xlim(testing_input[0], testing_input[-1])
+_ = plt.legend()
\ No newline at end of file
diff --git a/docs/source/examples.rst b/docs/source/examples.rst
index 2616ef02539e27094053f2f29b4947a164ecda9e..9f07b14c2282c362f790a502a2d57f1178040ada 100644
--- a/docs/source/examples.rst
+++ b/docs/source/examples.rst
@@ -3,4 +3,5 @@ Example Gallery
 
 .. toctree::
 
-   Emulator - RobustGaSP <../auto_examples/emulator/robustgasp/index>
\ No newline at end of file
+   Emulator - RobustGaSP <../auto_examples/emulator/robustgasp/index>
+   
\ No newline at end of file