[EXP] Reinitialization

examples/plot_estimate_gmm.py

@@ -20,12 +20,12 @@
 n_samples = 300
 n_features = 2
 X = np.ndarray((n_samples, n_features))
-X[:n_samples / 3, :] = random_state.multivariate_normal(
-    [0.0, 1.0], [[0.5, -1.0], [-1.0, 5.0]], size=(n_samples / 3,))
-X[n_samples / 3:-n_samples / 3, :] = random_state.multivariate_normal(
-    [-2.0, -2.0], [[3.0, 1.0], [1.0, 1.0]], size=(n_samples / 3,))
-X[-n_samples / 3:, :] = random_state.multivariate_normal(
-    [3.0, 1.0], [[3.0, -1.0], [-1.0, 1.0]], size=(n_samples / 3,))
+X[:n_samples // 3, :] = random_state.multivariate_normal(
+    [0.0, 1.0], [[0.5, -1.0], [-1.0, 5.0]], size=(n_samples // 3,))
+X[n_samples // 3:-n_samples // 3, :] = random_state.multivariate_normal(
+    [-2.0, -2.0], [[3.0, 1.0], [1.0, 1.0]], size=(n_samples // 3,))
+X[-n_samples // 3:, :] = random_state.multivariate_normal(
+    [3.0, 1.0], [[3.0, -1.0], [-1.0, 1.0]], size=(n_samples // 3,))
 
 gmm = GMM(n_components=3, random_state=random_state)
 gmm.from_samples(X)

examples/plot_initializations.py

@@ -0,0 +1,86 @@
+"""
+==================================
+Compare Initializations Strategies
+==================================
+
+Expectation Maximization for Gaussian Mixture Models does not have a unique
+solution. The result depends on the initialization. It is particularly
+important to either normalize the training data or set the covariances
+accordingly. In addition, k-means++ initialization helps to find a good
+initial distribution of means.
+"""
+print(__doc__)
+
+import numpy as np
+import matplotlib.pyplot as plt
+from gmr.utils import check_random_state
+from gmr import GMM, plot_error_ellipses, kmeansplusplus_initialization, covariance_initialization
+
+
+random_state = check_random_state(0)
+
+n_samples = 300
+n_features = 2
+X = np.ndarray((n_samples, n_features))
+X[:n_samples // 3, :] = random_state.multivariate_normal(
+    [0.0, 1.0], [[0.5, -1.0], [-1.0, 5.0]], size=(n_samples // 3,))
+X[n_samples // 3:-n_samples // 3, :] = random_state.multivariate_normal(
+    [-2.0, -2.0], [[3.0, 1.0], [1.0, 1.0]], size=(n_samples // 3,))
+X[-n_samples // 3:, :] = random_state.multivariate_normal(
+    [3.0, 1.0], [[3.0, -1.0], [-1.0, 1.0]], size=(n_samples // 3,))
+
+# artificial scaling, makes standard implementation fail
+# either the initial covariances have to be adjusted or we have
+# to normalize the dataset
+X[:, 1] *= 100.0
+
+plt.figure(figsize=(10, 10))
+
+n_components = 3
+initial_covs = np.empty((n_components, n_features, n_features))
+initial_covs[:] = np.eye(n_features)
+gmm = GMM(n_components=n_components,
+          priors=np.ones(n_components, dtype=np.float) / n_components,
+          means=np.zeros((n_components, n_features)),
+          covariances=initial_covs, random_state=random_state)
+
+plt.subplot(2, 2, 1)
+plt.title("Default initialization")
+plt.xlim((-10, 10))
+plt.ylim((-1000, 1000))
+plot_error_ellipses(plt.gca(), gmm, colors=["r", "g", "b"], alpha=0.15)
+plt.scatter(X[:, 0], X[:, 1])
+
+gmm.from_samples(X)
+
+plt.subplot(2, 2, 2)
+plt.title("Trained Gaussian Mixture Model")
+plt.xlim((-10, 10))
+plt.ylim((-1000, 1000))
+plot_error_ellipses(plt.gca(), gmm, colors=["r", "g", "b"], alpha=0.15)
+plt.scatter(X[:, 0], X[:, 1])
+
+initial_means = kmeansplusplus_initialization(X, n_components, random_state)
+initial_covs = covariance_initialization(X, n_components)
+gmm = GMM(n_components=n_components,
+          priors=np.ones(n_components, dtype=np.float) / n_components,
+          means=np.copy(initial_means),
+          covariances=initial_covs, random_state=random_state)
+
+plt.subplot(2, 2, 3)
+plt.title("k-means++ and inital covariance scaling")
+plt.xlim((-10, 10))
+plt.ylim((-1000, 1000))
+plot_error_ellipses(plt.gca(), gmm, colors=["r", "g", "b"], alpha=0.15)
+plt.scatter(X[:, 0], X[:, 1])
+
+gmm.from_samples(X)
+
+plt.subplot(2, 2, 4)
+plt.title("Trained Gaussian Mixture Model")
+plt.xlim((-10, 10))
+plt.ylim((-1000, 1000))
+plot_error_ellipses(plt.gca(), gmm, colors=["r", "g", "b"], alpha=0.15)
+plt.scatter(X[:, 0], X[:, 1])
+
+plt.show()

examples/plot_regression.py

@@ -20,7 +20,7 @@
 n_samples = 10
 X = np.ndarray((n_samples, 2))
 X[:, 0] = np.linspace(0, 2 * np.pi, n_samples)
-X[:, 1] = 1 - 3 * X[:, 0] + random_state.randn(n_samples)
+X[:, 1] = 100 * (1 - 3 * X[:, 0] + random_state.randn(n_samples))
 
 mvn = MVN(random_state=0)
 mvn.from_samples(X)
@@ -35,7 +35,7 @@
 plt.scatter(X[:, 0], X[:, 1])
 y = mean.ravel()
 s = covariance.ravel()
-plt.fill_between(X_test, y - s, y + s, alpha=0.2)
+plt.fill_between(X_test, y - np.sqrt(s), y + np.sqrt(s), alpha=0.2)
 plt.plot(X_test, y, lw=2)
 
 n_samples = 100

examples/plot_trajectories.py

@@ -0,0 +1,152 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.patches import Ellipse
+from itertools import cycle
+from gmr import GMM, plot_error_ellipses, kmeansplusplus_initialization, covariance_initialization, plot_error_ellipse
+from gmr.utils import check_random_state
+
+
+def make_demonstrations(n_demonstrations, n_steps, sigma=0.25, mu=0.5,
+                        start=np.zeros(2), goal=np.ones(2), random_state=None):
+    """Generates demonstration that can be used to test imitation learning.
+
+    Parameters
+    ----------
+    n_demonstrations : int
+        Number of noisy demonstrations
+
+    n_steps : int
+        Number of time steps
+
+    sigma : float, optional (default: 0.25)
+        Standard deviation of noisy component
+
+    mu : float, optional (default: 0.5)
+        Mean of noisy component
+
+    start : array, shape (2,), optional (default: 0s)
+        Initial position
+
+    goal : array, shape (2,), optional (default: 1s)
+        Final position
+
+    random_state : int
+        Seed for random number generator
+
+    Returns
+    -------
+    X : array, shape (n_task_dims, n_steps, n_demonstrations)
+        Noisy demonstrated trajectories
+
+    ground_truth : array, shape (n_task_dims, n_steps)
+        Original trajectory
+    """
+    random_state = np.random.RandomState(random_state)
+
+    X = np.empty((2, n_steps, n_demonstrations))
+
+    # Generate ground-truth for plotting
+    ground_truth = np.empty((2, n_steps))
+    T = np.linspace(-0, 1, n_steps)
+    ground_truth[0] = T
+    ground_truth[1] = (T / 20 + 1 / (sigma * np.sqrt(2 * np.pi)) *
+                       np.exp(-0.5 * ((T - mu) / sigma) ** 2))
+
+    # Generate trajectories
+    for i in range(n_demonstrations):
+        noisy_sigma = sigma * random_state.normal(1.0, 0.1)
+        noisy_mu = mu * random_state.normal(1.0, 0.1)
+        X[0, :, i] = T
+        X[1, :, i] = T + (1 / (noisy_sigma * np.sqrt(2 * np.pi)) *
+                          np.exp(-0.5 * ((T - noisy_mu) /
+                                         noisy_sigma) ** 2))
+
+    # Spatial alignment
+    current_start = ground_truth[:, 0]
+    current_goal = ground_truth[:, -1]
+    current_amplitude = current_goal - current_start
+    amplitude = goal - start
+    ground_truth = ((ground_truth.T - current_start) * amplitude /
+                    current_amplitude + start).T
+
+    for demo_idx in range(n_demonstrations):
+        current_start = X[:, 0, demo_idx]
+        current_goal = X[:, -1, demo_idx]
+        current_amplitude = current_goal - current_start
+        X[:, :, demo_idx] = ((X[:, :, demo_idx].T - current_start) *
+                             amplitude / current_amplitude + start).T
+
+    return X, ground_truth
+
+
+plot_covariances = False
+X, _ = make_demonstrations(
+    n_demonstrations=500, n_steps=20, goal=np.array([1., 2.]),
+    random_state=0)
+X = X.transpose(2, 1, 0)
+n_demonstrations, n_steps, n_task_dims = X.shape
+X_train = np.empty((n_demonstrations, n_steps, n_task_dims + 1))
+X_train[:, :, 1:] = X
+t = np.linspace(0, 1, n_steps)
+X_train[:, :, 0] = t
+X_train = X_train.reshape(n_demonstrations * n_steps, n_task_dims + 1)
+
+random_state = check_random_state(0)
+n_components = 10
+initial_means = kmeansplusplus_initialization(X_train, n_components, random_state)
+initial_covs = covariance_initialization(X_train, n_components)
+from sklearn.mixture import BayesianGaussianMixture
+bgmm = BayesianGaussianMixture(n_components=n_components, max_iter=500).fit(X_train)
+gmm = GMM(
+    n_components=n_components,
+    priors=bgmm.weights_,
+    means=bgmm.means_,
+    #means=np.copy(initial_means),
+    covariances=bgmm.covariances_,
+    #covariances=initial_covs,
+    verbose=2,
+    random_state=random_state)
+#gmm.from_samples(
+#    X_train, n_iter=100,
+#    reinit_means=True, min_eff_sample=n_task_dims, max_eff_sample=0.8)
+
+plt.figure()
+plt.subplot(111)
+
+plt.plot(X[:, :, 0].T, X[:, :, 1].T, c="k", alpha=0.1)
+
+means_over_time = []
+y_stds = []
+for step in t:
+    conditional_gmm = gmm.condition([0], np.array([step]))
+    conditional_mvn = conditional_gmm.to_mvn()
+    means_over_time.append(conditional_mvn.mean)
+    y_stds.append(conditional_mvn.covariance[1, 1])
+    samples = conditional_gmm.sample(100)
+    plt.scatter(samples[:, 0], samples[:, 1], s=1)
+means_over_time = np.array(means_over_time)
+y_stds = np.array(y_stds)
+plt.plot(means_over_time[:, 0], means_over_time[:, 1], c="r", lw=2)
+plt.fill_between(
+    means_over_time[:, 0],
+    means_over_time[:, 1] - 1.96 * y_stds,
+    means_over_time[:, 1] + 1.96 * y_stds,
+    color="r", alpha=0.5)
+
+if plot_covariances:
+    colors = cycle(["r", "g", "b"])
+    for factor in np.linspace(0.5, 4.0, 8):
+        new_gmm = GMM(
+            n_components=len(gmm.means), priors=gmm.priors,
+            means=gmm.means[:, 1:], covariances=gmm.covariances[:, 1:, 1:],
+            random_state=gmm.random_state)
+        for mean, (angle, width, height) in new_gmm.to_ellipses(factor):
+            ell = Ellipse(xy=mean, width=width, height=height,
+                            angle=np.degrees(angle))
+            ell.set_alpha(0.2)
+            ell.set_color(next(colors))
+            plt.gca().add_artist(ell)
+
+plt.xlabel("$x_1$")
+plt.ylabel("$x_2$")
+plt.show()

gmr/__init__.py

@@ -12,4 +12,5 @@
 __all__ = ['gmm', 'mvn', 'utils']
 
 from .mvn import MVN, plot_error_ellipse
-from .gmm import GMM, plot_error_ellipses
+from .gmm import (GMM, plot_error_ellipses, kmeansplusplus_initialization,
+                  covariance_initialization)