diff --git a/examples/models/gmm.py b/examples/models/gmm.py
deleted file mode 100644
index f46485e..0000000
--- a/examples/models/gmm.py
+++ /dev/null
@@ -1,84 +0,0 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-"""Provide some examples using GMMs.
-"""
-
-import numpy as np
-import matplotlib.pyplot as plt
-from sklearn.mixture import GaussianMixture
-
-from pyrobolearn.models.gmm import Gaussian, GMM, plot_gmm, plot_gmm_sklearn
-
-
-# create manually a GMM
-dim, num_components = 2, 5
-gmm = GMM(gaussians=[Gaussian(mean=np.random.uniform(-1., 1., size=dim),
-                              covariance=0.1*np.identity(dim)) for _ in range(num_components)])
-gmm_sklearn = GaussianMixture(n_components=num_components)
-
-
-# plot initial GMM
-plot_gmm(gmm, title='Initial GMM')
-plt.show()
-
-
-# create data: Generate random sample following a sine curve
-# Ref: https://scikit-learn.org/stable/auto_examples/mixture/plot_gmm_sin.html#sphx-glr-auto-examples-mixture-\
-# plot-gmm-sin-py
-n_samples = 100
-np.random.seed(0)
-X = np.zeros((n_samples, 2))
-step = 4. * np.pi / n_samples
-
-for i in range(X.shape[0]):
-    x = i * step - 6.
-    X[i, 0] = x + np.random.normal(0, 0.1)
-    X[i, 1] = 3. * (np.sin(x) + np.random.normal(0, .2))
-
-xlim, ylim = [-8, 8], [-8, 8]
-
-# plot data
-plt.title('Training data')
-plt.scatter(X[:, 0], X[:, 1])
-plt.show()
-
-
-# init GMM
-init_method = 'k-means'  # 'random', 'k-means', 'uniform', 'sklearn', 'curvature'
-gmm.init(X, method=init_method)
-fig, ax = plt.subplots(1, 1)
-plot_gmm(gmm, X=X, ax=ax, title='GMM after ' + init_method.capitalize(), xlim=xlim, ylim=ylim)
-plt.show()
-
-
-# fit a GMM using EM
-result = gmm.fit(X, init=None)
-gmm_sklearn.fit(X)
-
-# plot EM optimization
-plt.plot(result['losses'])
-plt.title('EM per iteration')
-plt.show()
-
-# plot trained GMM
-fig, ax = plt.subplots(1, 2)
-plot_gmm(gmm, X=X, label=True, ax=ax[0], title='Our Trained GMM', option=1, xlim=xlim, ylim=ylim)
-plot_gmm_sklearn(gmm_sklearn, X, label=True, ax=ax[1], title="Sklearn's Trained GMM", xlim=xlim, ylim=ylim)
-plt.show()
-
-# GMR: condition on the input variable and plot
-means, std_devs = [], []
-time_linspace = np.linspace(-6, 6, 100)
-for t in time_linspace:
-    g = gmm.condition(t, idx_out=1, idx_in=0).approximate_by_single_gaussian()
-    means.append(g.mean[0])
-    std_devs.append(np.sqrt(g.covariance[0, 0]))
-
-means, std_devs = np.asarray(means), np.asarray(std_devs)
-
-plt.plot(time_linspace, means)
-plt.fill_between(time_linspace, means - 2 * std_devs, means + 2 * std_devs, facecolor='green', alpha=0.3)
-plt.fill_between(time_linspace, means - std_devs, means + std_devs, facecolor='green', alpha=0.5)
-plt.title('GMR')
-plt.scatter(X[:, 0], X[:, 1])
-plt.show()
diff --git a/examples/models/gmr.py b/examples/models/gmr.py
index fccd1af..8a5f4f2 100644
--- a/examples/models/gmr.py
+++ b/examples/models/gmr.py
@@ -1,58 +1,59 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
-"""Provide some examples using GMR.
-
-See also the `gmm.py` example beforehand. In this example, we delve a bit deeper into GMR using 2D letters as training
-data.
+"""Provide some examples using GMM/GMR.
 """
 
 import numpy as np
-from scipy.io import loadmat
 import matplotlib.pyplot as plt
+from sklearn.mixture import GaussianMixture
 
-from pyrobolearn.models.gmm import Gaussian, GMM, plot_gmm, plot_gmr
+from pyrobolearn.models.gmm import Gaussian, GMM, plot_gmm, plot_gmm_sklearn
 
 
-# load the training data
-G = loadmat('../../data/2Dletters/G.mat')  # dict
-demos = G['demos']  # shape (1,N)
-n_demos = demos.shape[1]
-dim = demos[0, 0][0, 0][0].shape[0]
-length = demos[0, 0][0, 0][0].shape[1]
+# create manually a GMM
+dim, num_components = 2, 5
+gmm = GMM(gaussians=[Gaussian(mean=np.random.uniform(-1., 1., size=dim),
+                              covariance=0.1*np.identity(dim)) for _ in range(num_components)])
+gmm_sklearn = GaussianMixture(n_components=num_components)
 
-# plot the training data (x,y)
-X = []
-xlim, ylim = [-10, 10], [-10, 10]
-plt.xlim(xlim)
-plt.ylim(ylim)
-for i in range(0, n_demos, 2):
-    demo = demos[0, i][0, 0][0]  # shape (2, 200)
-    plt.plot(demo[0], demo[1])
-    X.append(demo.T)
+
+# plot initial GMM
+plot_gmm(gmm, title='Initial GMM')
 plt.show()
 
-# reshape training data (add time in addition to (x,y), thus we now have (t,x,y))
-time_linspace = np.linspace(0, 2., length)
-times = np.asarray([time_linspace for _ in range(len(X))]).reshape(-1, 1)
-X = np.vstack(X)  # shape (N*200, 2)
-X = np.hstack((times, X))  # shape (N*200, 3)
-print(X.shape)
 
-# create GMM
-dim, num_components = X.shape[1], 7
-gmm = GMM(gaussians=[Gaussian(mean=np.concatenate((np.random.uniform(0, 2., size=1),
-                                                   np.random.uniform(-8., 8., size=dim-1))),
-                              covariance=0.1*np.identity(dim)) for _ in range(num_components)])
+# create data: Generate random sample following a sine curve
+# Ref: https://scikit-learn.org/stable/auto_examples/mixture/plot_gmm_sin.html#sphx-glr-auto-examples-mixture-\
+# plot-gmm-sin-py
+n_samples = 100
+np.random.seed(0)
+X = np.zeros((n_samples, 2))
+step = 4. * np.pi / n_samples
+
+for i in range(X.shape[0]):
+    x = i * step - 6.
+    X[i, 0] = x + np.random.normal(0, 0.1)
+    X[i, 1] = 3. * (np.sin(x) + np.random.normal(0, .2))
+
+xlim, ylim = [-8, 8], [-8, 8]
+
+# plot data
+plt.title('Training data')
+plt.scatter(X[:, 0], X[:, 1])
+plt.show()
+
 
 # init GMM
 init_method = 'k-means'  # 'random', 'k-means', 'uniform', 'sklearn', 'curvature'
 gmm.init(X, method=init_method)
 fig, ax = plt.subplots(1, 1)
-plot_gmm(gmm, dims=[1, 2], X=X, ax=ax, title='GMM after ' + init_method.capitalize(), xlim=xlim, ylim=ylim)
+plot_gmm(gmm, X=X, ax=ax, title='GMM after ' + init_method.capitalize(), xlim=xlim, ylim=ylim)
 plt.show()
 
-# fit a GMM on it
-result = gmm.fit(X, init=None, num_iters=200)
+
+# fit a GMM using EM
+result = gmm.fit(X, init=None)
+gmm_sklearn.fit(X)
 
 # plot EM optimization
 plt.plot(result['losses'])
@@ -60,16 +61,24 @@ plt.title('EM per iteration')
 plt.show()
 
 # plot trained GMM
-fig, ax = plt.subplots(1, 1)
-plot_gmm(gmm, dims=[1, 2], X=X, label=True, ax=ax, title='Our Trained GMM', option=1, xlim=xlim, ylim=ylim)
+fig, ax = plt.subplots(1, 2)
+plot_gmm(gmm, X=X, label=True, ax=ax[0], title='Our Trained GMM', option=1, xlim=xlim, ylim=ylim)
+plot_gmm_sklearn(gmm_sklearn, X, label=True, ax=ax[1], title="Sklearn's Trained GMM", xlim=xlim, ylim=ylim)
 plt.show()
 
 # GMR: condition on the input variable and plot
-gaussians = []
+means, std_devs = [], []
+time_linspace = np.linspace(-6, 6, 100)
 for t in time_linspace:
-    g = gmm.condition(t, idx_out=[1, 2], idx_in=0).approximate_by_single_gaussian()
-    gaussians.append(g)
+    g = gmm.condition(t, idx_out=1, idx_in=0).approximate_by_single_gaussian()
+    means.append(g.mean[0])
+    std_devs.append(np.sqrt(g.covariance[0, 0]))
 
-# plot figures for GMR
-plot_gmr(time_linspace, gaussians=gaussians, xlim=xlim, ylim=ylim)
+means, std_devs = np.asarray(means), np.asarray(std_devs)
+
+plt.plot(time_linspace, means)
+plt.fill_between(time_linspace, means - 2 * std_devs, means + 2 * std_devs, facecolor='green', alpha=0.3)
+plt.fill_between(time_linspace, means - std_devs, means + std_devs, facecolor='green', alpha=0.5)
+plt.title('GMR')
+plt.scatter(X[:, 0], X[:, 1])
 plt.show()
diff --git a/examples/models/gmr_letter.py b/examples/models/gmr_letter.py
new file mode 100644
index 0000000..97fb687
--- /dev/null
+++ b/examples/models/gmr_letter.py
@@ -0,0 +1,76 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""Provide some examples using GMR.
+
+See also the `gmr.py` example beforehand. In this example, we delve a bit deeper into GMR using 2D letters as training
+data.
+"""
+
+import numpy as np
+from scipy.io import loadmat
+import matplotlib.pyplot as plt
+
+from pyrobolearn.models.gmm import Gaussian, GMM, plot_gmm, plot_gmr
+
+
+# load the training data
+G = loadmat('../../data/2Dletters/G.mat')  # dict
+demos = G['demos']  # shape (1,N)
+n_demos = demos.shape[1]
+dim = demos[0, 0][0, 0][0].shape[0]
+length = demos[0, 0][0, 0][0].shape[1]
+
+# plot the training data (x,y)
+X = []
+xlim, ylim = [-10, 10], [-10, 10]
+plt.title("Training Data")
+plt.xlim(xlim)
+plt.ylim(ylim)
+for i in range(0, n_demos, 2):
+    demo = demos[0, i][0, 0][0]  # shape (2, 200)
+    plt.plot(demo[0], demo[1])
+    X.append(demo.T)
+plt.show()
+
+# reshape training data (add time in addition to (x,y), thus we now have (t,x,y))
+time_linspace = np.linspace(0, 2., length)
+times = np.asarray([time_linspace for _ in range(len(X))]).reshape(-1, 1)
+X = np.vstack(X)  # shape (N*200, 2)
+X = np.hstack((times, X))  # shape (N*200, 3)
+print(X.shape)
+
+# create GMM
+dim, num_components = X.shape[1], 7
+gmm = GMM(gaussians=[Gaussian(mean=np.concatenate((np.random.uniform(0, 2., size=1),
+                                                   np.random.uniform(-8., 8., size=dim-1))),
+                              covariance=0.1*np.identity(dim)) for _ in range(num_components)])
+
+# init GMM
+init_method = 'k-means'  # 'random', 'k-means', 'uniform', 'sklearn', 'curvature'
+gmm.init(X, method=init_method)
+fig, ax = plt.subplots(1, 1)
+plot_gmm(gmm, dims=[1, 2], X=X, ax=ax, title='GMM after ' + init_method.capitalize(), xlim=xlim, ylim=ylim)
+plt.show()
+
+# fit a GMM on it
+result = gmm.fit(X, init=None, num_iters=200)
+
+# plot EM optimization
+plt.plot(result['losses'])
+plt.title('EM per iteration')
+plt.show()
+
+# plot trained GMM
+fig, ax = plt.subplots(1, 1)
+plot_gmm(gmm, dims=[1, 2], X=X, label=True, ax=ax, title='Our Trained GMM', option=1, xlim=xlim, ylim=ylim)
+plt.show()
+
+# GMR: condition on the input variable and plot
+gaussians = []
+for t in time_linspace:
+    g = gmm.condition(t, idx_out=[1, 2], idx_in=0).approximate_by_single_gaussian()
+    gaussians.append(g)
+
+# plot figures for GMR
+plot_gmr(time_linspace, gaussians=gaussians, xlim=xlim, ylim=ylim)
+plt.show()
diff --git a/examples/models/kmp_letter.py b/examples/models/kmp_letter.py
new file mode 100644
index 0000000..e933c39
--- /dev/null
+++ b/examples/models/kmp_letter.py
@@ -0,0 +1,82 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""Provide some examples using KMP.
+
+See also the `kmp.py` example beforehand. In this example, we delve a bit deeper into KMP using 2D letters as training
+data.
+"""
+
+import numpy as np
+from scipy.io import loadmat
+import matplotlib.pyplot as plt
+
+from pyrobolearn.models.gmm import plot_gmr, plot_gmm
+from pyrobolearn.models.kmp import KMP, RBF
+
+
+# KMP parameters (play with them)
+mean_reg = 1.  # 0.01, 0.1, 1.
+covariance_reg = 100.  # 0.1
+lengthscale = 1./6
+
+# load the training data
+G = loadmat('../../data/2Dletters/G.mat')  # dict
+demos = G['demos']  # shape (1,N)
+n_demos = demos.shape[1]
+dim = demos[0, 0][0, 0][0].shape[0]
+length = demos[0, 0][0, 0][0].shape[1]
+
+# plot the training data (x,y)
+X = []
+xlim, ylim = [-10, 10], [-10, 10]
+plt.title("Training Data")
+plt.xlim(xlim)
+plt.ylim(ylim)
+for i in range(0, n_demos, 2):
+    demo = demos[0, i][0, 0][0]  # shape (2, 200)
+    plt.plot(demo[0], demo[1])
+    X.append(demo.T)
+plt.show()
+
+# reshape training data (add time in addition to (x,y), thus we now have (t,x,y))
+time_linspace = np.linspace(0, 2., length)  # shape (200,)
+times = np.asarray([time_linspace for _ in range(len(X))])  # shape (N,200)
+X = np.asarray(X)  # shape (N, 200, 2)
+X = np.dstack((times, X))  # shape (N, 200, 3)
+print(X.shape)
+
+# create KMP
+print("Creating the KMP model")
+kernel = RBF(lengthscale=lengthscale)
+kmp = KMP(kernel_fct=kernel)
+
+# fit a KMP on the data
+print("Training the KMP...")
+kmp.fit(X=X[:, :, [0]], Y=X[:, :, 1:], gmm_num_components=7, mean_reg=mean_reg, covariance_reg=covariance_reg,
+        gmm_num_iters=200, database_size_limit=200, verbose=True)
+print("Finished the training")
+
+# plot underlying GMM
+gmm = kmp.reference_probability_distribution
+plot_gmm(gmm, dims=[1, 2], X=X.reshape(-1, 3), label=True, title='Underlying trained GMM', option=1, xlim=xlim,
+         ylim=ylim)
+plt.show()
+
+# predict with GMR
+gaussians = []
+for t in time_linspace:
+    g = gmm.condition(t, idx_out=[1, 2], idx_in=0).approximate_by_single_gaussian()
+    gaussians.append(g)
+
+# plot figures for GMR
+plot_gmr(time_linspace, gaussians=gaussians, xlim=xlim, ylim=ylim, suptitle='GMR')
+
+# predict with the KMP
+gaussians = []
+for t in time_linspace:
+    g = kmp.predict_proba(t, return_gaussian=True)
+    gaussians.append(g)
+
+# plot figures for KMP
+plot_gmr(time_linspace, gaussians=gaussians, xlim=xlim, ylim=ylim, suptitle='KMP')
+plt.show()
diff --git a/pyrobolearn/models/gmm/gmm.py b/pyrobolearn/models/gmm/gmm.py
index 5579f8d..90f2ab1 100755
--- a/pyrobolearn/models/gmm/gmm.py
+++ b/pyrobolearn/models/gmm/gmm.py
@@ -2146,6 +2146,7 @@ def plot_gmr(time_linspace, means=None, std_devs=None, covariances=None, gaussia
     plt.suptitle(suptitle)
 
     # plot t-x and t-y
+    limits = [xlim, ylim]
     for i in range(len(ylabels)):
         mean = means[:, i]
         std_dev = std_devs[:, i]
@@ -2154,6 +2155,8 @@ def plot_gmr(time_linspace, means=None, std_devs=None, covariances=None, gaussia
         axes[i].fill_between(time_linspace, mean - std_dev, mean + std_dev, facecolor='green', alpha=0.5)
         axes[i].set_xlabel('t')
         axes[i].set_ylabel(ylabels[i])
+        if i < len(limits):
+            axes[i].set_ylim(limits[i])
         # axes[i].scatter(X[:, 0], X[:, i+1])
 
     axes[-1].plot(means[:, 0], means[:, 1])
diff --git a/pyrobolearn/models/kmp/kmp.py b/pyrobolearn/models/kmp/kmp.py
index af05754..630218a 100644
--- a/pyrobolearn/models/kmp/kmp.py
+++ b/pyrobolearn/models/kmp/kmp.py
@@ -127,6 +127,9 @@ class KMP(object):
         else:
             self._database = database
 
+        # reference probability distribution (usually GMM)
+        self._ref_prob = None
+
         # mean
         self.mu = None  # mean used for the mean prediction
 
@@ -185,8 +188,8 @@ class KMP(object):
         if not isinstance(value, (int, float)):
             raise TypeError("Expecting the prior regularization term for the mean to be a scalar (int, float), but "
                             "got instead: {}".format(type(value)))
-        if value <= 0.:
-            raise ValueError("The prior regularization term needs to be strictly bigger than 0.")
+        if value < 0.:
+            raise ValueError("The prior regularization term for the mean needs to be bigger or equal to 0.")
         self._l1 = value
 
     # aliases
@@ -206,7 +209,8 @@ class KMP(object):
             raise TypeError("Expecting the prior regularization term for the covariance to be a scalar (int, float), "
                             "but got instead: {}".format(type(value)))
         if value <= 0.:
-            raise ValueError("The prior regularization term needs to be strictly bigger than 0.")
+            raise ValueError("The prior regularization term for the covariance needs to be strictly bigger than 0. "
+                             "The reason is that when computing the covariance prediction, we divide by that term.")
         self._l2 = value
 
     # aliases
@@ -214,6 +218,12 @@ class KMP(object):
     covariance_regularization = lambda2
     cov_reg = lambda2
 
+    @property
+    def reference_probability_distribution(self):
+        """Return the underlying reference probability distribution. Currently, this returns the underlying trained
+        GMM."""
+        return self._ref_prob
+
     ##################
     # Static Methods #
     ##################
@@ -310,6 +320,7 @@ class KMP(object):
         Returns:
             list[(np.ndarray, Gaussian)]: database which is a list of tuples where each one contains an input data
                 array and the corresponding predicted output Gaussian (by GMR)
+            GMM: reference probability distribution computed from the given data.
         """
         # TODO: replace gmm by joint generative model
 
@@ -384,10 +395,10 @@ class KMP(object):
                     for x in database]
 
         if verbose:
-            print("Database created...")
+            print("Database created with size: {}".format(len(database)))
 
         # return constructed reference database
-        return database
+        return database, gmm
 
     @staticmethod
     def get_reference_database(x, means=None, covariances=None, gaussians=None):
@@ -571,13 +582,16 @@ class KMP(object):
         # TODO: replace gmm by joint generative model
 
         # create reference database
-        self._database = self.create_reference_database(X, Y, gmm=gmm, gmm_num_components=gmm_num_components,
-                                                        distance=distance, database_threshold=database_threshold,
-                                                        database_size_limit=database_size_limit,
-                                                        sample_from_gmm=sample_from_gmm, gmm_init=gmm_init,
-                                                        gmm_reg=gmm_reg, gmm_num_iters=gmm_num_iters,
-                                                        gmm_convergence_threshold=gmm_convergence_threshold,
-                                                        seed=seed, verbose=verbose, block=block)
+        self._database, gmm = self.create_reference_database(X, Y, gmm=gmm, gmm_num_components=gmm_num_components,
+                                                             distance=distance, database_threshold=database_threshold,
+                                                             database_size_limit=database_size_limit,
+                                                             sample_from_gmm=sample_from_gmm, gmm_init=gmm_init,
+                                                             gmm_reg=gmm_reg, gmm_num_iters=gmm_num_iters,
+                                                             gmm_convergence_threshold=gmm_convergence_threshold,
+                                                             seed=seed, verbose=verbose, block=block)
+
+        # save reference probability distribution
+        self._ref_prob = gmm
 
         # compute kernel inverse from database
         K, K_inv1, K_inv2 = self.learn_from_database(self._database, mean_reg=mean_reg, covariance_reg=covariance_reg,