Merge branch 'devel' of github.com:SheffieldML/GPy into devel

GPy/examples/regression.py

@@ -1,7 +1,6 @@
 # Copyright (c) 2012, GPy authors (see AUTHORS.txt).
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
-
 """
 Gaussian Processes regression examples
 """
@@ -9,88 +8,107 @@ import pylab as pb
 import numpy as np
 import GPy
 
-def coregionalization_toy2(max_iters=100):
+def olympic_marathon_men(optimize=True, plot=True):
+    """Run a standard Gaussian process regression on the Olympic marathon data."""
+    data = GPy.util.datasets.olympic_marathon_men()
+
+    # create simple GP Model
+    m = GPy.models.GPRegression(data['X'], data['Y'])
+
+    # set the lengthscale to be something sensible (defaults to 1)
+    m['rbf_lengthscale'] = 10
+
+    if optimize:
+        m.optimize('bfgs', max_iters=200)
+    if plot:
+        m.plot(plot_limits=(1850, 2050))
+
+    return m
+
+def coregionalization_toy2(optimize=True, plot=True):
     """
     A simple demonstration of coregionalization on two sinusoidal functions.
     """
+    #build a design matrix with a column of integers indicating the output
     X1 = np.random.rand(50, 1) * 8
     X2 = np.random.rand(30, 1) * 5
     index = np.vstack((np.zeros_like(X1), np.ones_like(X2)))
     X = np.hstack((np.vstack((X1, X2)), index))
+
+    #build a suitable set of observed variables
     Y1 = np.sin(X1) + np.random.randn(*X1.shape) * 0.05
     Y2 = np.sin(X2) + np.random.randn(*X2.shape) * 0.05 + 2.
     Y = np.vstack((Y1, Y2))
 
+    #build the kernel
     k1 = GPy.kern.rbf(1) + GPy.kern.bias(1)
     k2 = GPy.kern.coregionalize(2,1)
-    k = k1**k2 #k = k1.prod(k2,tensor=True)
+    k = k1**k2
     m = GPy.models.GPRegression(X, Y, kernel=k)
     m.constrain_fixed('.*rbf_var', 1.)
-    # m.constrain_positive('.*kappa')
-    m.optimize('sim', messages=1, max_iters=max_iters)
 
-    pb.figure()
-    Xtest1 = np.hstack((np.linspace(0, 9, 100)[:, None], np.zeros((100, 1))))
-    Xtest2 = np.hstack((np.linspace(0, 9, 100)[:, None], np.ones((100, 1))))
-    mean, var, low, up = m.predict(Xtest1)
-    GPy.util.plot.gpplot(Xtest1[:, 0], mean, low, up)
-    mean, var, low, up = m.predict(Xtest2)
-    GPy.util.plot.gpplot(Xtest2[:, 0], mean, low, up)
-    pb.plot(X1[:, 0], Y1[:, 0], 'rx', mew=2)
-    pb.plot(X2[:, 0], Y2[:, 0], 'gx', mew=2)
+    if optimize:
+        m.optimize('bfgs', max_iters=100)
+
+    if plot:
+        m.plot(fixed_inputs=[(1,0)])
+        m.plot(fixed_inputs=[(1,1)], ax=pb.gca())
+
     return m
 
-def coregionalization_toy(max_iters=100):
-    """
-    A simple demonstration of coregionalization on two sinusoidal functions.
-    """
-    X1 = np.random.rand(50, 1) * 8
-    X2 = np.random.rand(30, 1) * 5
-    X = np.vstack((X1, X2))
-    Y1 = np.sin(X1) + np.random.randn(*X1.shape) * 0.05
-    Y2 = -np.sin(X2) + np.random.randn(*X2.shape) * 0.05
-    Y = np.vstack((Y1, Y2))
+#FIXME: Needs recovering once likelihoods are consolidated
+#def coregionalization_toy(optimize=True, plot=True):
+#    """
+#    A simple demonstration of coregionalization on two sinusoidal functions.
+#    """
+#    X1 = np.random.rand(50, 1) * 8
+#    X2 = np.random.rand(30, 1) * 5
+#    X = np.vstack((X1, X2))
+#    Y1 = np.sin(X1) + np.random.randn(*X1.shape) * 0.05
+#    Y2 = -np.sin(X2) + np.random.randn(*X2.shape) * 0.05
+#    Y = np.vstack((Y1, Y2))
+#
+#    k1 = GPy.kern.rbf(1)
+#    m = GPy.models.GPMultioutputRegression(X_list=[X1,X2],Y_list=[Y1,Y2],kernel_list=[k1])
+#    m.constrain_fixed('.*rbf_var', 1.)
+#    m.optimize(max_iters=100)
+#
+#    fig, axes = pb.subplots(2,1)
+#    m.plot(fixed_inputs=[(1,0)],ax=axes[0])
+#    m.plot(fixed_inputs=[(1,1)],ax=axes[1])
+#    axes[0].set_title('Output 0')
+#    axes[1].set_title('Output 1')
+#    return m
 
-    k1 = GPy.kern.rbf(1)
-    m = GPy.models.GPMultioutputRegression(X_list=[X1,X2],Y_list=[Y1,Y2],kernel_list=[k1])
-    m.constrain_fixed('.*rbf_var', 1.)
-    m.optimize(max_iters=max_iters)
-
-    fig, axes = pb.subplots(2,1)
-    m.plot(fixed_inputs=[(1,0)],ax=axes[0])
-    m.plot(fixed_inputs=[(1,1)],ax=axes[1])
-    axes[0].set_title('Output 0')
-    axes[1].set_title('Output 1')
-    return m
-
-def coregionalization_sparse(max_iters=100):
+def coregionalization_sparse(optimize=True, plot=True):
     """
     A simple demonstration of coregionalization on two sinusoidal functions using sparse approximations.
     """
-    X1 = np.random.rand(500, 1) * 8
-    X2 = np.random.rand(300, 1) * 5
-    index = np.vstack((np.zeros_like(X1), np.ones_like(X2)))
-    X = np.hstack((np.vstack((X1, X2)), index))
-    Y1 = np.sin(X1) + np.random.randn(*X1.shape) * 0.05
-    Y2 = -np.sin(X2) + np.random.randn(*X2.shape) * 0.05
-    Y = np.vstack((Y1, Y2))
+    #fetch the data from the non sparse examples
+    m = coregionalization_toy2(optimize=False, plot=False)
+    X, Y = m.X, m.likelihood.Y
 
-    k1 = GPy.kern.rbf(1)
+    #construct a model
+    m = GPy.models.SparseGPRegression(X,Y)
+    m.constrain_fixed('iip_\d+_1') # don't optimize the inducing input indexes
 
-    m = GPy.models.SparseGPMultioutputRegression(X_list=[X1,X2],Y_list=[Y1,Y2],kernel_list=[k1],num_inducing=5)
-    m.constrain_fixed('.*rbf_var',1.)
-    #m.optimize(messages=1)
-    m.optimize_restarts(5, robust=True, messages=1, max_iters=max_iters, optimizer='bfgs')
+    if optimize:
+        m.optimize('bfgs', max_iters=100, messages=1)
+
+    if plot:
+        m.plot(fixed_inputs=[(1,0)])
+        m.plot(fixed_inputs=[(1,1)], ax=pb.gca())
 
-    fig, axes = pb.subplots(2,1)
-    m.plot_single_output(output=0,ax=axes[0],plot_limits=(-1,9))
-    m.plot_single_output(output=1,ax=axes[1],plot_limits=(-1,9))
-    axes[0].set_title('Output 0')
-    axes[1].set_title('Output 1')
     return m
 
-def epomeo_gpx(max_iters=100):
-    """Perform Gaussian process regression on the latitude and longitude data from the Mount Epomeo runs. Requires gpxpy to be installed on your system to load in the data."""
+
+
+def epomeo_gpx(optimize=True, plot=True):
+    """
+    Perform Gaussian process regression on the latitude and longitude data
+    from the Mount Epomeo runs. Requires gpxpy to be installed on your system
+    to load in the data.
+    """
     data = GPy.util.datasets.epomeo_gpx()
     num_data_list = []
     for Xpart in data['X']:
@@ -119,14 +137,17 @@ def epomeo_gpx(max_iters=100):
     m.constrain_fixed('.*rbf_var', 1.)
     m.constrain_fixed('iip')
     m.constrain_bounded('noise_variance', 1e-3, 1e-1)
-#     m.optimize_restarts(5, robust=True, messages=1, max_iters=max_iters, optimizer='bfgs')
     m.optimize(max_iters=max_iters,messages=True)
 
     return m
 
 
 def multiple_optima(gene_number=937, resolution=80, model_restarts=10, seed=10000, max_iters=300):
-    """Show an example of a multimodal error surface for Gaussian process regression. Gene 939 has bimodal behaviour where the noisy mode is higher."""
+    """
+    Show an example of a multimodal error surface for Gaussian process
+    regression. Gene 939 has bimodal behaviour where the noisy mode is
+    higher.
+    """
 
     # Contour over a range of length scales and signal/noise ratios.
     length_scales = np.linspace(0.1, 60., resolution)
@@ -175,12 +196,15 @@ def multiple_optima(gene_number=937, resolution=80, model_restarts=10, seed=1000
     return m # (models, lls)
 
 def _contour_data(data, length_scales, log_SNRs, kernel_call=GPy.kern.rbf):
-    """Evaluate the GP objective function for a given data set for a range of signal to noise ratios and a range of lengthscales.
+    """
+    Evaluate the GP objective function for a given data set for a range of
+    signal to noise ratios and a range of lengthscales.
 
     :data_set: A data set from the utils.datasets director.
     :length_scales: a list of length scales to explore for the contour plot.
     :log_SNRs: a list of base 10 logarithm signal to noise ratios to explore for the contour plot.
-    :kernel: a kernel to use for the 'signal' portion of the data."""
+    :kernel: a kernel to use for the 'signal' portion of the data.
+    """
 
     lls = []
     total_var = np.var(data['Y'])
@@ -203,79 +227,58 @@ def _contour_data(data, length_scales, log_SNRs, kernel_call=GPy.kern.rbf):
     return np.array(lls)
 
 
-def olympic_100m_men(max_iters=100, kernel=None):
+def olympic_100m_men(optimize=True, plot=True):
     """Run a standard Gaussian process regression on the Rogers and Girolami olympics data."""
     data = GPy.util.datasets.olympic_100m_men()
 
     # create simple GP Model
-    m = GPy.models.GPRegression(data['X'], data['Y'], kernel)
+    m = GPy.models.GPRegression(data['X'], data['Y'])
 
     # set the lengthscale to be something sensible (defaults to 1)
-    if kernel==None:
-        m['rbf_lengthscale'] = 10
+    m['rbf_lengthscale'] = 10
 
-    # optimize
-    m.optimize(max_iters=max_iters)
+    if optimize:
+        m.optimize('bfgs', max_iters=200)
 
-    # plot
-    m.plot(plot_limits=(1850, 2050))
-    print(m)
+    if plot:
+        m.plot(plot_limits=(1850, 2050))
     return m
 
-def olympic_marathon_men(max_iters=100, kernel=None):
-    """Run a standard Gaussian process regression on the Olympic marathon data."""
-    data = GPy.util.datasets.olympic_marathon_men()
-
-    # create simple GP Model
-    m = GPy.models.GPRegression(data['X'], data['Y'], kernel)
-
-    # set the lengthscale to be something sensible (defaults to 1)
-    if kernel==None:
-        m['rbf_lengthscale'] = 10
-
-    # optimize
-    m.optimize(max_iters=max_iters)
-
-    # plot
-    m.plot(plot_limits=(1850, 2050))
-    print(m)
-    return m
-
-def toy_rbf_1d(optimizer='tnc', max_nb_eval_optim=100):
+def toy_rbf_1d(optimize=True, plot=True):
     """Run a simple demonstration of a standard Gaussian process fitting it to data sampled from an RBF covariance."""
     data = GPy.util.datasets.toy_rbf_1d()
 
     # create simple GP Model
     m = GPy.models.GPRegression(data['X'], data['Y'])
 
-    # optimize
-    m.optimize(optimizer, max_f_eval=max_nb_eval_optim)
-    # plot
-    m.plot()
-    print(m)
+    if optimize:
+        m.optimize('bfgs')
+    if plot:
+        m.plot()
+
     return m
 
-def toy_rbf_1d_50(max_iters=100):
+def toy_rbf_1d_50(optimize=True, plot=True):
     """Run a simple demonstration of a standard Gaussian process fitting it to data sampled from an RBF covariance."""
     data = GPy.util.datasets.toy_rbf_1d_50()
 
     # create simple GP Model
     m = GPy.models.GPRegression(data['X'], data['Y'])
 
-    # optimize
-    m.optimize(max_iters=max_iters)
+    if optimize:
+        m.optimize('bfgs')
+    if plot:
+        m.plot()
 
-    # plot
-    m.plot()
-    print(m)
     return m
 
-def toy_poisson_rbf_1d(optimizer='bfgs', max_nb_eval_optim=100):
+
+def toy_poisson_rbf_1d(optimize=True, plot=True):
     """Run a simple demonstration of a standard Gaussian process fitting it to data sampled from an RBF covariance."""
     x_len = 400
     X = np.linspace(0, 10, x_len)[:, None]
     f_true = np.random.multivariate_normal(np.zeros(x_len), GPy.kern.rbf(1).K(X))
-    Y = np.array([np.random.poisson(np.exp(f)) for f in f_true])[:,None]
+    Y = np.array([np.random.poisson(np.exp(f)) for f in f_true]).reshape(x_len,1)
 
     noise_model = GPy.likelihoods.poisson()
     likelihood = GPy.likelihoods.EP(Y,noise_model)
@@ -283,14 +286,14 @@ def toy_poisson_rbf_1d(optimizer='bfgs', max_nb_eval_optim=100):
     # create simple GP Model
     m = GPy.models.GPRegression(X, Y, likelihood=likelihood)
 
-    # optimize
-    m.optimize(optimizer, max_f_eval=max_nb_eval_optim)
-    # plot
-    m.plot()
-    print(m)
+    if optimize:
+        m.optimize('bfgs')
+    if plot:
+        m.plot()
+
     return m
 
-def toy_poisson_rbf_1d_laplace(optimizer='bfgs', max_nb_eval_optim=100):
+def toy_poisson_rbf_1d_laplace(optimize=True, plot=True):
     """Run a simple demonstration of a standard Gaussian process fitting it to data sampled from an RBF covariance."""
     x_len = 30
     X = np.linspace(0, 10, x_len)[:, None]
@@ -303,13 +306,13 @@ def toy_poisson_rbf_1d_laplace(optimizer='bfgs', max_nb_eval_optim=100):
     # create simple GP Model
     m = GPy.models.GPRegression(X, Y, likelihood=likelihood)
 
-    # optimize
-    m.optimize(optimizer, max_f_eval=max_nb_eval_optim)
-    # plot
-    m.plot()
-    # plot the real underlying rate function
-    pb.plot(X, np.exp(f_true), '--k', linewidth=2)
-    print(m)
+    if optimize:
+        m.optimize(optimizer, max_f_eval=max_nb_eval_optim)
+    if plot:
+        m.plot()
+        # plot the real underlying rate function
+        pb.plot(X, np.exp(f_true), '--k', linewidth=2)
+
     return m
 
 
@@ -459,7 +462,7 @@ def sparse_GP_regression_2D(num_samples=400, num_inducing=50, max_iters=100):
     print(m)
     return m
 
-def uncertain_inputs_sparse_regression(max_iters=100):
+def uncertain_inputs_sparse_regression(optimize=True, plot=True):
     """Run a 1D example of a sparse GP regression with uncertain inputs."""
     fig, axes = pb.subplots(1, 2, figsize=(12, 5))
 

GPy/models.py

@@ -1,23 +1,31 @@
 '''
-Created on 14 Nov 2013
+GPy Models
+==========
 
-@author: maxz
+Implementations for common models used in GP regression and classification.
+The different models can be viewed in :mod:`GPy.models_modules`, which holds
+detailed explanations for the different models.
+
+:warning: This module is a convienince module for endusers to use. For developers 
+see :mod:`GPy.models_modules`, which holds the implementions for each model. 
 '''
 
-from _models.bayesian_gplvm import BayesianGPLVM
-from _models.gp_regression import GPRegression
-from _models.gp_classification import GPClassification#; _gp_classification = gp_classification ; del gp_classification 
-from _models.sparse_gp_regression import SparseGPRegression#; _sparse_gp_regression = sparse_gp_regression ; del sparse_gp_regression 
-from _models.svigp_regression import SVIGPRegression#; _svigp_regression = svigp_regression ; del svigp_regression 
-from _models.sparse_gp_classification import SparseGPClassification#; _sparse_gp_classification = sparse_gp_classification ; del sparse_gp_classification 
-from _models.fitc_classification import FITCClassification#; _fitc_classification = fitc_classification ; del fitc_classification 
-from _models.gplvm import GPLVM#; _gplvm = gplvm ; del gplvm 
-from _models.bcgplvm import BCGPLVM#; _bcgplvm = bcgplvm; del bcgplvm
-from _models.sparse_gplvm import SparseGPLVM#; _sparse_gplvm = sparse_gplvm ; del sparse_gplvm 
-from _models.warped_gp import WarpedGP#; _warped_gp = warped_gp ; del warped_gp 
-from _models.bayesian_gplvm import BayesianGPLVM#; _bayesian_gplvm = bayesian_gplvm ; del bayesian_gplvm 
-from _models.mrd import MRD#; _mrd = mrd; del mrd 
-from _models.gradient_checker import GradientChecker#; _gradient_checker = gradient_checker ; del gradient_checker 
-from _models.gp_multioutput_regression import GPMultioutputRegression#; _gp_multioutput_regression = gp_multioutput_regression ; del gp_multioutput_regression 
-from _models.sparse_gp_multioutput_regression import SparseGPMultioutputRegression#; _sparse_gp_multioutput_regression = sparse_gp_multioutput_regression ; del sparse_gp_multioutput_regression 
-from _models.gradient_checker import GradientChecker
+__updated__ = '2013-11-28'
+
+from models_modules.bayesian_gplvm import BayesianGPLVM
+from models_modules.gp_regression import GPRegression
+from models_modules.gp_classification import GPClassification#; _gp_classification = gp_classification ; del gp_classification 
+from models_modules.sparse_gp_regression import SparseGPRegression#; _sparse_gp_regression = sparse_gp_regression ; del sparse_gp_regression 
+from models_modules.svigp_regression import SVIGPRegression#; _svigp_regression = svigp_regression ; del svigp_regression 
+from models_modules.sparse_gp_classification import SparseGPClassification#; _sparse_gp_classification = sparse_gp_classification ; del sparse_gp_classification 
+from models_modules.fitc_classification import FITCClassification#; _fitc_classification = fitc_classification ; del fitc_classification 
+from models_modules.gplvm import GPLVM#; _gplvm = gplvm ; del gplvm 
+from models_modules.bcgplvm import BCGPLVM#; _bcgplvm = bcgplvm; del bcgplvm
+from models_modules.sparse_gplvm import SparseGPLVM#; _sparse_gplvm = sparse_gplvm ; del sparse_gplvm 
+from models_modules.warped_gp import WarpedGP#; _warped_gp = warped_gp ; del warped_gp 
+from models_modules.bayesian_gplvm import BayesianGPLVM#; _bayesian_gplvm = bayesian_gplvm ; del bayesian_gplvm 
+from models_modules.mrd import MRD#; _mrd = mrd; del mrd 
+from models_modules.gradient_checker import GradientChecker#; _gradient_checker = gradient_checker ; del gradient_checker 
+from models_modules.gp_multioutput_regression import GPMultioutputRegression#; _gp_multioutput_regression = gp_multioutput_regression ; del gp_multioutput_regression 
+from models_modules.sparse_gp_multioutput_regression import SparseGPMultioutputRegression#; _sparse_gp_multioutput_regression = sparse_gp_multioutput_regression ; del sparse_gp_multioutput_regression 
+from models_modules.gradient_checker import GradientChecker

GPy/models_modules/bayesian_gplvm.py

@@ -12,6 +12,7 @@ from GPy.util import plot_latent, linalg
 from .gplvm import GPLVM
 from GPy.util.plot_latent import most_significant_input_dimensions
 from matplotlib import pyplot
+from GPy.core.model import Model
 
 class BayesianGPLVM(SparseGP, GPLVM):
     """
@@ -285,6 +286,57 @@ class BayesianGPLVM(SparseGP, GPLVM):
         self.init = state.pop()
         SparseGP.setstate(self, state)
 
+class BayesianGPLVMWithMissingData(Model):
+    """
+    Bayesian Gaussian Process Latent Variable Model with missing data support.
+    NOTE: Missing data is assumed to be missing at random!
+    
+    This extension comes with a large memory and computing time deficiency.
+    Use only if fraction of missing data at random is higher than 60%.
+    Otherwise, try filtering data before using this extension.
+    
+    Y can hold missing data as given by `missing`, standard is :class:`~numpy.nan`.
+    
+    If likelihood is given for Y, this likelihood will be discarded, but the parameters
+    of the likelihood will be taken. Also every effort of creating the same likelihood
+    will be done.
+     
+    :param likelihood_or_Y: observed data (np.ndarray) or GPy.likelihood
+    :type likelihood_or_Y: :class:`~numpy.ndarray` | :class:`~GPy.likelihoods.likelihood.likelihood` instance
+    :param int input_dim: latent dimensionality
+    :param init: initialisation method for the latent space
+    :type init: 'PCA' | 'random'
+    """
+    def __init__(self, likelihood_or_Y, input_dim, X=None, X_variance=None, init='PCA', num_inducing=10,
+                 Z=None, kernel=None, missing=np.nan, **kwargs):
+        if type(likelihood_or_Y) is np.ndarray:
+            likelihood = Gaussian(likelihood_or_Y)
+        else:
+            likelihood = likelihood_or_Y
+
+        if X == None:
+            X = self.initialise_latent(init, input_dim, likelihood.Y)
+        self.init = init
+
+        if X_variance is None:
+            X_variance = np.clip((np.ones_like(X) * 0.5) + .01 * np.random.randn(*X.shape), 0.001, 1)
+
+        if Z is None:
+            Z = np.random.permutation(X.copy())[:num_inducing]
+        assert Z.shape[1] == X.shape[1]
+
+        if kernel is None:
+            kernel = kern.rbf(input_dim) # + kern.white(input_dim)
+
+        SparseGP.__init__(self, X, likelihood, kernel, Z=Z, X_variance=X_variance, **kwargs)
+        self.ensure_default_constraints()
+
+    def _get_param_names(self):
+        X_names = sum([['X_%i_%i' % (n, q) for q in range(self.input_dim)] for n in range(self.num_data)], [])
+        S_names = sum([['X_variance_%i_%i' % (n, q) for q in range(self.input_dim)] for n in range(self.num_data)], [])
+        return (X_names + S_names + SparseGP._get_param_names(self))
+
+    pass
 
 def latent_cost_and_grad(mu_S, kern, Z, dL_dpsi0, dL_dpsi1, dL_dpsi2):
     """