Added optimize and plot for classification, non_gaussian and stochastic examples

GPy/examples/classification.py

@@ -6,12 +6,11 @@
 Gaussian Processes classification
 """
 import pylab as pb
-import numpy as np
 import GPy
 
 default_seed = 10000
 
-def oil(num_inducing=50, max_iters=100, kernel=None):
+def oil(num_inducing=50, max_iters=100, kernel=None, optimize=True, plot=True):
     """
     Run a Gaussian process classification on the three phase oil data. The demonstration calls the basic GP classification model and uses EP to approximate the likelihood.
 
@@ -25,7 +24,7 @@ def oil(num_inducing=50, max_iters=100, kernel=None):
     Ytest[Ytest.flatten()==-1] = 0
 
     # Create GP model
-    m = GPy.models.SparseGPClassification(X, Y,kernel=kernel,num_inducing=num_inducing)
+    m = GPy.models.SparseGPClassification(X, Y, kernel=kernel, num_inducing=num_inducing)
 
     # Contrain all parameters to be positive
     m.tie_params('.*len')
@@ -33,15 +32,16 @@ def oil(num_inducing=50, max_iters=100, kernel=None):
     m.update_likelihood_approximation()
 
     # Optimize
-    m.optimize(max_iters=max_iters)
+    if optimize:
+        m.optimize(max_iters=max_iters)
     print(m)
 
     #Test
     probs = m.predict(Xtest)[0]
-    GPy.util.classification.conf_matrix(probs,Ytest)
+    GPy.util.classification.conf_matrix(probs, Ytest)
     return m
 
-def toy_linear_1d_classification(seed=default_seed):
+def toy_linear_1d_classification(seed=default_seed, optimize=True, plot=True):
     """
     Simple 1D classification example using EP approximation
 
@@ -58,21 +58,23 @@ def toy_linear_1d_classification(seed=default_seed):
     m = GPy.models.GPClassification(data['X'], Y)
 
     # Optimize
-    #m.update_likelihood_approximation()
+    if optimize:
-    # Parameters optimization:
+        #m.update_likelihood_approximation()
-    #m.optimize()
+        # Parameters optimization:
-    #m.update_likelihood_approximation()
+        #m.optimize()
-    m.pseudo_EM()
+        #m.update_likelihood_approximation()
+        m.pseudo_EM()
 
     # Plot
-    fig, axes = pb.subplots(2,1)
+    if plot:
-    m.plot_f(ax=axes[0])
+        fig, axes = pb.subplots(2, 1)
-    m.plot(ax=axes[1])
+        m.plot_f(ax=axes[0])
-    print(m)
+        m.plot(ax=axes[1])
 
+    print m
     return m
 
-def toy_linear_1d_classification_laplace(seed=default_seed):
+def toy_linear_1d_classification_laplace(seed=default_seed, optimize=True, plot=True):
     """
     Simple 1D classification example using Laplace approximation
 
@@ -90,24 +92,25 @@ def toy_linear_1d_classification_laplace(seed=default_seed):
 
     # Model definition
     m = GPy.models.GPClassification(data['X'], Y, likelihood=laplace_likelihood)
-
     print m
+
     # Optimize
-    #m.update_likelihood_approximation()
+    if optimize:
-    # Parameters optimization:
+        #m.update_likelihood_approximation()
-    m.optimize('bfgs', messages=1)
+        # Parameters optimization:
-    #m.pseudo_EM()
+        m.optimize('bfgs', messages=1)
+        #m.pseudo_EM()
 
     # Plot
-    fig, axes = pb.subplots(2,1)
+    if plot:
-    m.plot_f(ax=axes[0])
+        fig, axes = pb.subplots(2, 1)
-    m.plot(ax=axes[1])
+        m.plot_f(ax=axes[0])
-    print(m)
+        m.plot(ax=axes[1])
 
+    print m
     return m
 
-
+def sparse_toy_linear_1d_classification(num_inducing=10, seed=default_seed, optimize=True, plot=True):
-def sparse_toy_linear_1d_classification(num_inducing=10,seed=default_seed):
     """
     Sparse 1D classification example
 
@@ -121,24 +124,26 @@ def sparse_toy_linear_1d_classification(num_inducing=10,seed=default_seed):
     Y[Y.flatten() == -1] = 0
 
     # Model definition
-    m = GPy.models.SparseGPClassification(data['X'], Y,num_inducing=num_inducing)
+    m = GPy.models.SparseGPClassification(data['X'], Y, num_inducing=num_inducing)
-    m['.*len']= 4.
+    m['.*len'] = 4.
 
     # Optimize
-    #m.update_likelihood_approximation()
+    if optimize:
-    # Parameters optimization:
+        #m.update_likelihood_approximation()
-    #m.optimize()
+        # Parameters optimization:
-    m.pseudo_EM()
+        #m.optimize()
+        m.pseudo_EM()
 
     # Plot
-    fig, axes = pb.subplots(2,1)
+    if plot:
-    m.plot_f(ax=axes[0])
+        fig, axes = pb.subplots(2, 1)
-    m.plot(ax=axes[1])
+        m.plot_f(ax=axes[0])
-    print(m)
+        m.plot(ax=axes[1])
 
+    print m
     return m
 
-def toy_heaviside(seed=default_seed):
+def toy_heaviside(seed=default_seed, optimize=True, plot=True):
     """
     Simple 1D classification example using a heavy side gp transformation
 
@@ -153,24 +158,26 @@ def toy_heaviside(seed=default_seed):
 
     # Model definition
     noise_model = GPy.likelihoods.bernoulli(GPy.likelihoods.noise_models.gp_transformations.Heaviside())
-    likelihood = GPy.likelihoods.EP(Y,noise_model)
+    likelihood = GPy.likelihoods.EP(Y, noise_model)
     m = GPy.models.GPClassification(data['X'], likelihood=likelihood)
 
     # Optimize
-    m.update_likelihood_approximation()
+    if optimize:
-    # Parameters optimization:
+        m.update_likelihood_approximation()
-    m.optimize()
+        # Parameters optimization:
-    #m.pseudo_EM()
+        m.optimize()
+        #m.pseudo_EM()
 
     # Plot
-    fig, axes = pb.subplots(2,1)
+    if plot:
-    m.plot_f(ax=axes[0])
+        fig, axes = pb.subplots(2, 1)
-    m.plot(ax=axes[1])
+        m.plot_f(ax=axes[0])
-    print(m)
+        m.plot(ax=axes[1])
 
+    print m
     return m
 
-def crescent_data(model_type='Full', num_inducing=10, seed=default_seed, kernel=None):
+def crescent_data(model_type='Full', num_inducing=10, seed=default_seed, kernel=None, optimize=True, plot=True):
     """
     Run a Gaussian process classification on the crescent data. The demonstration calls the basic GP classification model and uses EP to approximate the likelihood.
 
@@ -187,7 +194,7 @@ def crescent_data(model_type='Full', num_inducing=10, seed=default_seed, kernel=
     Y[Y.flatten()==-1] = 0
 
     if model_type == 'Full':
-        m = GPy.models.GPClassification(data['X'], Y,kernel=kernel)
+        m = GPy.models.GPClassification(data['X'], Y, kernel=kernel)
 
     elif model_type == 'DTC':
         m = GPy.models.SparseGPClassification(data['X'], Y, kernel=kernel, num_inducing=num_inducing)
@@ -197,8 +204,11 @@ def crescent_data(model_type='Full', num_inducing=10, seed=default_seed, kernel=
         m = GPy.models.FITCClassification(data['X'], Y, kernel=kernel, num_inducing=num_inducing)
         m['.*len'] = 3.
 
-    m.pseudo_EM()
+    if optimize:
-    print(m)
+        m.pseudo_EM()
-    m.plot()
 
+    if plot:
+        m.plot()
+
+    print m
     return m

GPy/examples/non_gaussian.py

@@ -114,7 +114,7 @@ def student_t_approx(optimize=True, plot=True):
 
     return m1, m2, m3, m4
 
-def boston_example():
+def boston_example(optimize=True, plot=True):
     import sklearn
     from sklearn.cross_validation import KFold
     optimizer='bfgs'
@@ -143,7 +143,6 @@ def boston_example():
         noise = 1e-1 #np.exp(-2)
         rbf_len = 0.5
         data_axis_plot = 4
-        plot = False
         kernelstu = GPy.kern.rbf(X.shape[1]) + GPy.kern.white(X.shape[1]) + GPy.kern.bias(X.shape[1])
         kernelgp = GPy.kern.rbf(X.shape[1]) + GPy.kern.white(X.shape[1]) + GPy.kern.bias(X.shape[1])
 
@@ -158,17 +157,13 @@ def boston_example():
         mgp['rbf_len'] = rbf_len
         mgp['noise'] = noise
         print mgp
-        mgp.optimize(optimizer=optimizer, messages=messages)
+        if optimize:
+            mgp.optimize(optimizer=optimizer, messages=messages)
         Y_test_pred = mgp.predict(X_test)
         score_folds[1, n] = rmse(Y_test, Y_test_pred[0])
         pred_density[1, n] = np.mean(mgp.log_predictive_density(X_test, Y_test))
         print mgp
         print pred_density
-        if plot:
-            plt.figure()
-            plt.scatter(X_test[:, data_axis_plot], Y_test_pred[0])
-            plt.scatter(X_test[:, data_axis_plot], Y_test, c='r', marker='x')
-            plt.title('GP gauss')
 
         print "Gaussian Laplace GP"
         N, D = Y_train.shape
@@ -181,20 +176,13 @@ def boston_example():
         mg['rbf_len'] = rbf_len
         mg['noise'] = noise
         print mg
-        try:
+        if optimize:
             mg.optimize(optimizer=optimizer, messages=messages)
-        except Exception:
-            print "Blew up"
         Y_test_pred = mg.predict(X_test)
         score_folds[2, n] = rmse(Y_test, Y_test_pred[0])
         pred_density[2, n] = np.mean(mg.log_predictive_density(X_test, Y_test))
         print pred_density
         print mg
-        if plot:
-            plt.figure()
-            plt.scatter(X_test[:, data_axis_plot], Y_test_pred[0])
-            plt.scatter(X_test[:, data_axis_plot], Y_test, c='r', marker='x')
-            plt.title('Lap gauss')
 
         for stu_num, df in enumerate(degrees_freedoms):
             #Student T
@@ -208,61 +196,71 @@ def boston_example():
             mstu_t['rbf_len'] = rbf_len
             mstu_t['t_noise'] = noise
             print mstu_t
-            try:
+            if optimize:
                 mstu_t.optimize(optimizer=optimizer, messages=messages)
-            except Exception:
-                print "Blew up"
             Y_test_pred = mstu_t.predict(X_test)
             score_folds[3+stu_num, n] = rmse(Y_test, Y_test_pred[0])
             pred_density[3+stu_num, n] = np.mean(mstu_t.log_predictive_density(X_test, Y_test))
             print pred_density
             print mstu_t
-            if plot:
+
-                plt.figure()
+    if plot:
-                plt.scatter(X_test[:, data_axis_plot], Y_test_pred[0])
+        plt.figure()
-                plt.scatter(X_test[:, data_axis_plot], Y_test, c='r', marker='x')
+        plt.scatter(X_test[:, data_axis_plot], Y_test_pred[0])
-                plt.title('Stu t {}df'.format(df))
+        plt.scatter(X_test[:, data_axis_plot], Y_test, c='r', marker='x')
+        plt.title('GP gauss')
+
+        plt.figure()
+        plt.scatter(X_test[:, data_axis_plot], Y_test_pred[0])
+        plt.scatter(X_test[:, data_axis_plot], Y_test, c='r', marker='x')
+        plt.title('Lap gauss')
+
+        plt.figure()
+        plt.scatter(X_test[:, data_axis_plot], Y_test_pred[0])
+        plt.scatter(X_test[:, data_axis_plot], Y_test, c='r', marker='x')
+        plt.title('Stu t {}df'.format(df))
 
     print "Average scores: {}".format(np.mean(score_folds, 1))
     print "Average pred density: {}".format(np.mean(pred_density, 1))
 
-    #Plotting
+    if plot:
-    stu_t_legends = ['Student T, df={}'.format(df) for df in degrees_freedoms]
+        #Plotting
-    legends = ['Baseline', 'Gaussian', 'Laplace Approx Gaussian'] + stu_t_legends
+        stu_t_legends = ['Student T, df={}'.format(df) for df in degrees_freedoms]
+        legends = ['Baseline', 'Gaussian', 'Laplace Approx Gaussian'] + stu_t_legends
 
-    #Plot boxplots for RMSE density
+        #Plot boxplots for RMSE density
-    fig = plt.figure()
+        fig = plt.figure()
-    ax=fig.add_subplot(111)
+        ax=fig.add_subplot(111)
-    plt.title('RMSE')
+        plt.title('RMSE')
-    bp = ax.boxplot(score_folds.T, notch=0, sym='+', vert=1, whis=1.5)
+        bp = ax.boxplot(score_folds.T, notch=0, sym='+', vert=1, whis=1.5)
-    plt.setp(bp['boxes'], color='black')
+        plt.setp(bp['boxes'], color='black')
-    plt.setp(bp['whiskers'], color='black')
+        plt.setp(bp['whiskers'], color='black')
-    plt.setp(bp['fliers'], color='red', marker='+')
+        plt.setp(bp['fliers'], color='red', marker='+')
-    xtickNames = plt.setp(ax, xticklabels=legends)
+        xtickNames = plt.setp(ax, xticklabels=legends)
-    plt.setp(xtickNames, rotation=45, fontsize=8)
+        plt.setp(xtickNames, rotation=45, fontsize=8)
-    ax.set_ylabel('RMSE')
+        ax.set_ylabel('RMSE')
-    ax.set_xlabel('Distribution')
+        ax.set_xlabel('Distribution')
-    #Make grid and put it below boxes
+        #Make grid and put it below boxes
-    ax.yaxis.grid(True, linestyle='-', which='major', color='lightgrey',
+        ax.yaxis.grid(True, linestyle='-', which='major', color='lightgrey',
-              alpha=0.5)
+                alpha=0.5)
-    ax.set_axisbelow(True)
+        ax.set_axisbelow(True)
 
-    #Plot boxplots for predictive density
+        #Plot boxplots for predictive density
-    fig = plt.figure()
+        fig = plt.figure()
-    ax=fig.add_subplot(111)
+        ax=fig.add_subplot(111)
-    plt.title('Predictive density')
+        plt.title('Predictive density')
-    bp = ax.boxplot(pred_density[1:,:].T, notch=0, sym='+', vert=1, whis=1.5)
+        bp = ax.boxplot(pred_density[1:,:].T, notch=0, sym='+', vert=1, whis=1.5)
-    plt.setp(bp['boxes'], color='black')
+        plt.setp(bp['boxes'], color='black')
-    plt.setp(bp['whiskers'], color='black')
+        plt.setp(bp['whiskers'], color='black')
-    plt.setp(bp['fliers'], color='red', marker='+')
+        plt.setp(bp['fliers'], color='red', marker='+')
-    xtickNames = plt.setp(ax, xticklabels=legends[1:])
+        xtickNames = plt.setp(ax, xticklabels=legends[1:])
-    plt.setp(xtickNames, rotation=45, fontsize=8)
+        plt.setp(xtickNames, rotation=45, fontsize=8)
-    ax.set_ylabel('Mean Log probability P(Y*|Y)')
+        ax.set_ylabel('Mean Log probability P(Y*|Y)')
-    ax.set_xlabel('Distribution')
+        ax.set_xlabel('Distribution')
-    #Make grid and put it below boxes
+        #Make grid and put it below boxes
-    ax.yaxis.grid(True, linestyle='-', which='major', color='lightgrey',
+        ax.yaxis.grid(True, linestyle='-', which='major', color='lightgrey',
-              alpha=0.5)
+                alpha=0.5)
-    ax.set_axisbelow(True)
+        ax.set_axisbelow(True)
     return mstu_t
 
 def precipitation_example():

GPy/examples/stochastic.py

@@ -5,7 +5,7 @@ import pylab as pb
 import numpy as np
 import GPy
 
-def toy_1d():
+def toy_1d(optimize=True, plot=True):
     N = 2000
     M = 20
 
@@ -20,15 +20,18 @@ def toy_1d():
 
     m.param_steplength = 1e-4
 
-    fig = pb.figure()
+    if plot:
-    ax = fig.add_subplot(111)
+        fig = pb.figure()
-    def cb():
+        ax = fig.add_subplot(111)
-        ax.cla()
+        def cb(foo):
-        m.plot(ax=ax,Z_height=-3)
+            ax.cla()
-        ax.set_ylim(-3,3)
+            m.plot(ax=ax,Z_height=-3)
-        fig.canvas.draw()
+            ax.set_ylim(-3,3)
+            fig.canvas.draw()
 
-    m.optimize(500, callback=cb, callback_interval=1)
+    if optimize:
+        m.optimize(500, callback=cb, callback_interval=1)
 
-    m.plot_traces()
+    if plot:
+        m.plot_traces()
     return m