mirror of
https://github.com/SheffieldML/GPy.git
synced 2026-05-10 04:22:38 +02:00
Huge merge correcting upstream master
This commit is contained in:
Alan Saul
commit
34932f8746
319 changed files with 26201 additions and 26660 deletions
|
|
@ -1,8 +1,7 @@
|
|||
# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
|
||||
# Copyright (c) 2012-2014, GPy authors (see AUTHORS.txt).
|
||||
# Licensed under the BSD 3-clause license (see LICENSE.txt)
|
||||
|
||||
import classification
|
||||
import regression
|
||||
import dimensionality_reduction
|
||||
import tutorials
|
||||
import stochastic
|
||||
import non_gaussian
|
||||
|
|
|
|||
|
|
@ -1,11 +1,10 @@
|
|||
# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
|
||||
# Copyright (c) 2012-2014, GPy authors (see AUTHORS.txt).
|
||||
# Licensed under the BSD 3-clause license (see LICENSE.txt)
|
||||
|
||||
|
||||
"""
|
||||
Gaussian Processes classification
|
||||
Gaussian Processes classification examples
|
||||
"""
|
||||
import pylab as pb
|
||||
import GPy
|
||||
|
||||
default_seed = 10000
|
||||
|
|
@ -15,7 +14,9 @@ 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.
|
||||
|
||||
"""
|
||||
data = GPy.util.datasets.oil()
|
||||
try:import pods
|
||||
except ImportError:print 'pods unavailable, see https://github.com/sods/ods for example datasets'
|
||||
data = pods.datasets.oil()
|
||||
X = data['X']
|
||||
Xtest = data['Xtest']
|
||||
Y = data['Y'][:, 0:1]
|
||||
|
|
@ -27,13 +28,13 @@ def oil(num_inducing=50, max_iters=100, kernel=None, optimize=True, plot=True):
|
|||
m = GPy.models.SparseGPClassification(X, Y, kernel=kernel, num_inducing=num_inducing)
|
||||
|
||||
# Contrain all parameters to be positive
|
||||
m.tie_params('.*len')
|
||||
#m.tie_params('.*len')
|
||||
m['.*len'] = 10.
|
||||
m.update_likelihood_approximation()
|
||||
|
||||
# Optimize
|
||||
if optimize:
|
||||
m.optimize(max_iters=max_iters)
|
||||
for _ in range(5):
|
||||
m.optimize(max_iters=int(max_iters/5))
|
||||
print(m)
|
||||
|
||||
#Test
|
||||
|
|
@ -50,7 +51,9 @@ def toy_linear_1d_classification(seed=default_seed, optimize=True, plot=True):
|
|||
|
||||
"""
|
||||
|
||||
data = GPy.util.datasets.toy_linear_1d_classification(seed=seed)
|
||||
try:import pods
|
||||
except ImportError:print 'pods unavailable, see https://github.com/sods/ods for example datasets'
|
||||
data = pods.datasets.toy_linear_1d_classification(seed=seed)
|
||||
Y = data['Y'][:, 0:1]
|
||||
Y[Y.flatten() == -1] = 0
|
||||
|
||||
|
|
@ -61,13 +64,14 @@ def toy_linear_1d_classification(seed=default_seed, optimize=True, plot=True):
|
|||
if optimize:
|
||||
#m.update_likelihood_approximation()
|
||||
# Parameters optimization:
|
||||
#m.optimize()
|
||||
m.optimize()
|
||||
#m.update_likelihood_approximation()
|
||||
m.pseudo_EM()
|
||||
#m.pseudo_EM()
|
||||
|
||||
# Plot
|
||||
if plot:
|
||||
fig, axes = pb.subplots(2, 1)
|
||||
from matplotlib import pyplot as plt
|
||||
fig, axes = plt.subplots(2, 1)
|
||||
m.plot_f(ax=axes[0])
|
||||
m.plot(ax=axes[1])
|
||||
|
||||
|
|
@ -83,27 +87,30 @@ def toy_linear_1d_classification_laplace(seed=default_seed, optimize=True, plot=
|
|||
|
||||
"""
|
||||
|
||||
data = GPy.util.datasets.toy_linear_1d_classification(seed=seed)
|
||||
try:import pods
|
||||
except ImportError:print 'pods unavailable, see https://github.com/sods/ods for example datasets'
|
||||
data = pods.datasets.toy_linear_1d_classification(seed=seed)
|
||||
Y = data['Y'][:, 0:1]
|
||||
Y[Y.flatten() == -1] = 0
|
||||
|
||||
bern_noise_model = GPy.likelihoods.bernoulli()
|
||||
laplace_likelihood = GPy.likelihoods.Laplace(Y.copy(), bern_noise_model)
|
||||
likelihood = GPy.likelihoods.Bernoulli()
|
||||
laplace_inf = GPy.inference.latent_function_inference.Laplace()
|
||||
kernel = GPy.kern.RBF(1)
|
||||
|
||||
# Model definition
|
||||
m = GPy.models.GPClassification(data['X'], Y, likelihood=laplace_likelihood)
|
||||
print m
|
||||
m = GPy.core.GP(data['X'], Y, kernel=kernel, likelihood=likelihood, inference_method=laplace_inf)
|
||||
|
||||
# Optimize
|
||||
if optimize:
|
||||
#m.update_likelihood_approximation()
|
||||
# Parameters optimization:
|
||||
m.optimize('bfgs', messages=1)
|
||||
#m.pseudo_EM()
|
||||
try:
|
||||
m.optimize('scg', messages=1)
|
||||
except Exception as e:
|
||||
return m
|
||||
|
||||
# Plot
|
||||
if plot:
|
||||
fig, axes = pb.subplots(2, 1)
|
||||
from matplotlib import pyplot as plt
|
||||
fig, axes = plt.subplots(2, 1)
|
||||
m.plot_f(ax=axes[0])
|
||||
m.plot(ax=axes[1])
|
||||
|
||||
|
|
@ -119,7 +126,9 @@ def sparse_toy_linear_1d_classification(num_inducing=10, seed=default_seed, opti
|
|||
|
||||
"""
|
||||
|
||||
data = GPy.util.datasets.toy_linear_1d_classification(seed=seed)
|
||||
try:import pods
|
||||
except ImportError:print 'pods unavailable, see https://github.com/sods/ods for example datasets'
|
||||
data = pods.datasets.toy_linear_1d_classification(seed=seed)
|
||||
Y = data['Y'][:, 0:1]
|
||||
Y[Y.flatten() == -1] = 0
|
||||
|
||||
|
|
@ -129,21 +138,19 @@ def sparse_toy_linear_1d_classification(num_inducing=10, seed=default_seed, opti
|
|||
|
||||
# Optimize
|
||||
if optimize:
|
||||
#m.update_likelihood_approximation()
|
||||
# Parameters optimization:
|
||||
#m.optimize()
|
||||
m.pseudo_EM()
|
||||
m.optimize()
|
||||
|
||||
# Plot
|
||||
if plot:
|
||||
fig, axes = pb.subplots(2, 1)
|
||||
from matplotlib import pyplot as plt
|
||||
fig, axes = plt.subplots(2, 1)
|
||||
m.plot_f(ax=axes[0])
|
||||
m.plot(ax=axes[1])
|
||||
|
||||
print m
|
||||
return m
|
||||
|
||||
def toy_heaviside(seed=default_seed, optimize=True, plot=True):
|
||||
def toy_heaviside(seed=default_seed, max_iters=100, optimize=True, plot=True):
|
||||
"""
|
||||
Simple 1D classification example using a heavy side gp transformation
|
||||
|
||||
|
|
@ -152,25 +159,30 @@ def toy_heaviside(seed=default_seed, optimize=True, plot=True):
|
|||
|
||||
"""
|
||||
|
||||
data = GPy.util.datasets.toy_linear_1d_classification(seed=seed)
|
||||
try:import pods
|
||||
except ImportError:print 'pods unavailable, see https://github.com/sods/ods for example datasets'
|
||||
data = pods.datasets.toy_linear_1d_classification(seed=seed)
|
||||
Y = data['Y'][:, 0:1]
|
||||
Y[Y.flatten() == -1] = 0
|
||||
|
||||
# Model definition
|
||||
noise_model = GPy.likelihoods.bernoulli(GPy.likelihoods.noise_models.gp_transformations.Heaviside())
|
||||
likelihood = GPy.likelihoods.EP(Y, noise_model)
|
||||
m = GPy.models.GPClassification(data['X'], likelihood=likelihood)
|
||||
kernel = GPy.kern.RBF(1)
|
||||
likelihood = GPy.likelihoods.Bernoulli(gp_link=GPy.likelihoods.link_functions.Heaviside())
|
||||
ep = GPy.inference.latent_function_inference.expectation_propagation.EP()
|
||||
m = GPy.core.GP(X=data['X'], Y=Y, kernel=kernel, likelihood=likelihood, inference_method=ep, name='gp_classification_heaviside')
|
||||
#m = GPy.models.GPClassification(data['X'], likelihood=likelihood)
|
||||
|
||||
# Optimize
|
||||
if optimize:
|
||||
m.update_likelihood_approximation()
|
||||
# Parameters optimization:
|
||||
m.optimize()
|
||||
#m.pseudo_EM()
|
||||
for _ in range(5):
|
||||
m.optimize(max_iters=int(max_iters/5))
|
||||
print m
|
||||
|
||||
# Plot
|
||||
if plot:
|
||||
fig, axes = pb.subplots(2, 1)
|
||||
from matplotlib import pyplot as plt
|
||||
fig, axes = plt.subplots(2, 1)
|
||||
m.plot_f(ax=axes[0])
|
||||
m.plot(ax=axes[1])
|
||||
|
||||
|
|
@ -189,7 +201,9 @@ def crescent_data(model_type='Full', num_inducing=10, seed=default_seed, kernel=
|
|||
:param kernel: kernel to use in the model
|
||||
:type kernel: a GPy kernel
|
||||
"""
|
||||
data = GPy.util.datasets.crescent_data(seed=seed)
|
||||
try:import pods
|
||||
except ImportError:print 'pods unavailable, see https://github.com/sods/ods for example datasets'
|
||||
data = pods.datasets.crescent_data(seed=seed)
|
||||
Y = data['Y']
|
||||
Y[Y.flatten()==-1] = 0
|
||||
|
||||
|
|
|
|||
89
GPy/examples/coreg_example.py
Normal file
89
GPy/examples/coreg_example.py
Normal file
|
|
@ -0,0 +1,89 @@
|
|||
# Copyright (c) 2012-2014, GPy authors (see AUTHORS.txt).
|
||||
# Licensed under the BSD 3-clause license (see LICENSE.txt)
|
||||
|
||||
import numpy as np
|
||||
try:
|
||||
import pylab as pb
|
||||
except:
|
||||
pass
|
||||
import GPy
|
||||
pb.ion()
|
||||
pb.close('all')
|
||||
|
||||
X1 = np.arange(3)[:,None]
|
||||
X2 = np.arange(4)[:,None]
|
||||
I1 = np.zeros_like(X1)
|
||||
I2 = np.ones_like(X2)
|
||||
|
||||
_X = np.vstack([ X1, X2 ])
|
||||
_I = np.vstack([ I1, I2 ])
|
||||
|
||||
X = np.hstack([ _X, _I ])
|
||||
|
||||
Y1 = np.sin(X1/8.)
|
||||
Y2 = np.cos(X2/8.)
|
||||
|
||||
Bias = GPy.kern.Bias(1,active_dims=[0])
|
||||
Coreg = GPy.kern.Coregionalize(1,2,active_dims=[1])
|
||||
K = Bias.prod(Coreg,name='X')
|
||||
|
||||
#K.coregion.W = 0
|
||||
#print K.coregion.W
|
||||
#print Bias.K(_X,_X)
|
||||
#print K.K(X,X)
|
||||
#pb.matshow(K.K(X,X))
|
||||
|
||||
Mlist = [GPy.kern.Matern32(1,lengthscale=20.,name="Mat")]
|
||||
kern = GPy.util.multioutput.LCM(input_dim=1,num_outputs=2,kernels_list=Mlist,name='H')
|
||||
kern.B.W = 0
|
||||
kern.B.kappa = 1.
|
||||
#kern.B.W.fix()
|
||||
#kern.B.kappa.fix()
|
||||
#m = GPy.models.GPCoregionalizedRegression(X_list=[X1,X2], Y_list=[Y1,Y2], kernel=kern)
|
||||
|
||||
|
||||
Z1 = np.array([1.5,2.5])[:,None]
|
||||
|
||||
m = GPy.models.SparseGPCoregionalizedRegression(X_list=[X1], Y_list=[Y1], Z_list = [Z1], kernel=kern)
|
||||
#m.optimize()
|
||||
m.checkgrad(verbose=1)
|
||||
|
||||
"""
|
||||
fig = pb.figure()
|
||||
ax0 = fig.add_subplot(211)
|
||||
ax1 = fig.add_subplot(212)
|
||||
slices = GPy.util.multioutput.get_slices([Y1,Y2])
|
||||
m.plot(fixed_inputs=[(1,0)],which_data_rows=slices[0],ax=ax0)
|
||||
#m.plot(fixed_inputs=[(1,1)],which_data_rows=slices[1],ax=ax1)
|
||||
"""
|
||||
|
||||
|
||||
|
||||
"""
|
||||
|
||||
X1 = 100 * np.random.rand(100)[:,None]
|
||||
X2 = 100 * np.random.rand(100)[:,None]
|
||||
#X1.sort()
|
||||
#X2.sort()
|
||||
|
||||
Y1 = np.sin(X1/10.) + np.random.rand(100)[:,None]
|
||||
Y2 = np.cos(X2/10.) + np.random.rand(100)[:,None]
|
||||
|
||||
|
||||
|
||||
|
||||
Mlist = [GPy.kern.Matern32(1,lengthscale=20.,name="Mat")]
|
||||
kern = GPy.util.multioutput.LCM(input_dim=1,num_outputs=12,kernels_list=Mlist,name='H')
|
||||
|
||||
|
||||
m = GPy.models.GPCoregionalizedRegression(X_list=[X1,X2], Y_list=[Y1,Y2], kernel=kern)
|
||||
m.optimize()
|
||||
|
||||
fig = pb.figure()
|
||||
ax0 = fig.add_subplot(211)
|
||||
ax1 = fig.add_subplot(212)
|
||||
slices = GPy.util.multioutput.get_slices([Y1,Y2])
|
||||
m.plot(fixed_inputs=[(1,0)],which_data_rows=slices[0],ax=ax0)
|
||||
m.plot(fixed_inputs=[(1,1)],which_data_rows=slices[1],ax=ax1)
|
||||
|
||||
"""
|
||||
|
|
@ -1,75 +1,80 @@
|
|||
# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
|
||||
# Copyright (c) 2012-2014, GPy authors (see AUTHORS.txt).
|
||||
# Licensed under the BSD 3-clause license (see LICENSE.txt)
|
||||
import numpy as _np
|
||||
default_seed = 123344
|
||||
|
||||
def bgplvm_test_model(seed=default_seed, optimize=False, verbose=1, plot=False):
|
||||
# default_seed = _np.random.seed(123344)
|
||||
|
||||
def bgplvm_test_model(optimize=False, verbose=1, plot=False, output_dim=200, nan=False):
|
||||
"""
|
||||
model for testing purposes. Samples from a GP with rbf kernel and learns
|
||||
the samples with a new kernel. Normally not for optimization, just model cheking
|
||||
"""
|
||||
from GPy.likelihoods.gaussian import Gaussian
|
||||
import GPy
|
||||
|
||||
num_inputs = 13
|
||||
num_inducing = 5
|
||||
if plot:
|
||||
output_dim = 1
|
||||
input_dim = 2
|
||||
input_dim = 3
|
||||
else:
|
||||
input_dim = 2
|
||||
output_dim = 25
|
||||
output_dim = output_dim
|
||||
|
||||
# generate GPLVM-like data
|
||||
X = _np.random.rand(num_inputs, input_dim)
|
||||
lengthscales = _np.random.rand(input_dim)
|
||||
k = (GPy.kern.rbf(input_dim, .5, lengthscales, ARD=True)
|
||||
+ GPy.kern.white(input_dim, 0.01))
|
||||
k = GPy.kern.RBF(input_dim, .5, lengthscales, ARD=True)
|
||||
K = k.K(X)
|
||||
Y = _np.random.multivariate_normal(_np.zeros(num_inputs), K, output_dim).T
|
||||
lik = Gaussian(Y, normalize=True)
|
||||
Y = _np.random.multivariate_normal(_np.zeros(num_inputs), K, (output_dim,)).T
|
||||
|
||||
k = GPy.kern.rbf_inv(input_dim, .5, _np.ones(input_dim) * 2., ARD=True) + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim)
|
||||
# k = GPy.kern.linear(input_dim) + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||
# k = GPy.kern.rbf(input_dim, ARD = False) + GPy.kern.white(input_dim, 0.00001)
|
||||
# k = GPy.kern.rbf(input_dim, .5, _np.ones(input_dim) * 2., ARD=True) + GPy.kern.rbf(input_dim, .3, _np.ones(input_dim) * .2, ARD=True)
|
||||
# k = GPy.kern.rbf(input_dim, .5, 2., ARD=0) + GPy.kern.rbf(input_dim, .3, .2, ARD=0)
|
||||
# k = GPy.kern.rbf(input_dim, .5, _np.ones(input_dim) * 2., ARD=True) + GPy.kern.linear(input_dim, _np.ones(input_dim) * .2, ARD=True)
|
||||
# k = GPy.kern.RBF_inv(input_dim, .5, _np.ones(input_dim) * 2., ARD=True) + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim)
|
||||
# k = GPy.kern.linear(input_dim)# + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||
# k = GPy.kern.RBF(input_dim, ARD = False) + GPy.kern.white(input_dim, 0.00001)
|
||||
# k = GPy.kern.RBF(input_dim, .5, _np.ones(input_dim) * 2., ARD=True) + GPy.kern.RBF(input_dim, .3, _np.ones(input_dim) * .2, ARD=True)
|
||||
# k = GPy.kern.RBF(input_dim, .5, 2., ARD=0) + GPy.kern.RBF(input_dim, .3, .2, ARD=0)
|
||||
# k = GPy.kern.RBF(input_dim, .5, _np.ones(input_dim) * 2., ARD=True) + GPy.kern.linear(input_dim, _np.ones(input_dim) * .2, ARD=True)
|
||||
|
||||
p = .3
|
||||
|
||||
m = GPy.models.BayesianGPLVM(Y, input_dim, kernel=k, num_inducing=num_inducing)
|
||||
|
||||
if nan:
|
||||
m.inference_method = GPy.inference.latent_function_inference.var_dtc.VarDTCMissingData()
|
||||
m.Y[_np.random.binomial(1, p, size=(Y.shape)).astype(bool)] = _np.nan
|
||||
m.parameters_changed()
|
||||
|
||||
m = GPy.models.BayesianGPLVM(lik, input_dim, kernel=k, num_inducing=num_inducing)
|
||||
#===========================================================================
|
||||
# randomly obstruct data with percentage p
|
||||
p = .8
|
||||
Y_obstruct = Y.copy()
|
||||
Y_obstruct[_np.random.uniform(size=(Y.shape)) < p] = _np.nan
|
||||
#===========================================================================
|
||||
m2 = GPy.models.BayesianGPLVMWithMissingData(Y_obstruct, input_dim, kernel=k, num_inducing=num_inducing)
|
||||
m.lengthscales = lengthscales
|
||||
# m2 = GPy.models.BayesianGPLVMWithMissingData(Y_obstruct, input_dim, kernel=k, num_inducing=num_inducing)
|
||||
# m.lengthscales = lengthscales
|
||||
|
||||
if plot:
|
||||
import matplotlib.pyplot as pb
|
||||
m.plot()
|
||||
pb.title('PCA initialisation')
|
||||
m2.plot()
|
||||
pb.title('PCA initialisation')
|
||||
# m2.plot()
|
||||
# pb.title('PCA initialisation')
|
||||
|
||||
if optimize:
|
||||
m.optimize('scg', messages=verbose)
|
||||
m2.optimize('scg', messages=verbose)
|
||||
# m2.optimize('scg', messages=verbose)
|
||||
if plot:
|
||||
m.plot()
|
||||
pb.title('After optimisation')
|
||||
m2.plot()
|
||||
pb.title('After optimisation')
|
||||
# m2.plot()
|
||||
# pb.title('After optimisation')
|
||||
|
||||
return m, m2
|
||||
return m
|
||||
|
||||
def gplvm_oil_100(optimize=True, verbose=1, plot=True):
|
||||
import GPy
|
||||
data = GPy.util.datasets.oil_100()
|
||||
import pods
|
||||
data = pods.datasets.oil_100()
|
||||
Y = data['X']
|
||||
# create simple GP model
|
||||
kernel = GPy.kern.rbf(6, ARD=True) + GPy.kern.bias(6)
|
||||
kernel = GPy.kern.RBF(6, ARD=True) + GPy.kern.Bias(6)
|
||||
m = GPy.models.GPLVM(Y, 6, kernel=kernel)
|
||||
m.data_labels = data['Y'].argmax(axis=1)
|
||||
if optimize: m.optimize('scg', messages=verbose)
|
||||
|
|
@ -78,13 +83,15 @@ def gplvm_oil_100(optimize=True, verbose=1, plot=True):
|
|||
|
||||
def sparse_gplvm_oil(optimize=True, verbose=0, plot=True, N=100, Q=6, num_inducing=15, max_iters=50):
|
||||
import GPy
|
||||
import pods
|
||||
|
||||
_np.random.seed(0)
|
||||
data = GPy.util.datasets.oil()
|
||||
data = pods.datasets.oil()
|
||||
Y = data['X'][:N]
|
||||
Y = Y - Y.mean(0)
|
||||
Y /= Y.std(0)
|
||||
# Create the model
|
||||
kernel = GPy.kern.rbf(Q, ARD=True) + GPy.kern.bias(Q)
|
||||
kernel = GPy.kern.RBF(Q, ARD=True) + GPy.kern.Bias(Q)
|
||||
m = GPy.models.SparseGPLVM(Y, Q, kernel=kernel, num_inducing=num_inducing)
|
||||
m.data_labels = data['Y'][:N].argmax(axis=1)
|
||||
|
||||
|
|
@ -94,9 +101,9 @@ def sparse_gplvm_oil(optimize=True, verbose=0, plot=True, N=100, Q=6, num_induci
|
|||
m.kern.plot_ARD()
|
||||
return m
|
||||
|
||||
def swiss_roll(optimize=True, verbose=1, plot=True, N=1000, num_inducing=15, Q=4, sigma=.2):
|
||||
def swiss_roll(optimize=True, verbose=1, plot=True, N=1000, num_inducing=25, Q=4, sigma=.2):
|
||||
import GPy
|
||||
from GPy.util.datasets import swiss_roll_generated
|
||||
from pods.datasets import swiss_roll_generated
|
||||
from GPy.models import BayesianGPLVM
|
||||
|
||||
data = swiss_roll_generated(num_samples=N, sigma=sigma)
|
||||
|
|
@ -134,93 +141,103 @@ def swiss_roll(optimize=True, verbose=1, plot=True, N=1000, num_inducing=15, Q=4
|
|||
(1 - var))) + .001
|
||||
Z = _np.random.permutation(X)[:num_inducing]
|
||||
|
||||
kernel = GPy.kern.rbf(Q, ARD=True) + GPy.kern.bias(Q, _np.exp(-2)) + GPy.kern.white(Q, _np.exp(-2))
|
||||
kernel = GPy.kern.RBF(Q, ARD=True) + GPy.kern.Bias(Q, _np.exp(-2)) + GPy.kern.White(Q, _np.exp(-2))
|
||||
|
||||
m = BayesianGPLVM(Y, Q, X=X, X_variance=S, num_inducing=num_inducing, Z=Z, kernel=kernel)
|
||||
m.data_colors = c
|
||||
m.data_t = t
|
||||
m['noise_variance'] = Y.var() / 100.
|
||||
|
||||
if optimize:
|
||||
m.optimize('scg', messages=verbose, max_iters=2e3)
|
||||
m.optimize('bfgs', messages=verbose, max_iters=2e3)
|
||||
|
||||
if plot:
|
||||
fig = plt.figure('fitted')
|
||||
ax = fig.add_subplot(111)
|
||||
s = m.input_sensitivity().argsort()[::-1][:2]
|
||||
ax.scatter(*m.X.T[s], c=c)
|
||||
ax.scatter(*m.X.mean.T[s], c=c)
|
||||
|
||||
return m
|
||||
|
||||
def bgplvm_oil(optimize=True, verbose=1, plot=True, N=200, Q=7, num_inducing=40, max_iters=1000, **k):
|
||||
import GPy
|
||||
from GPy.likelihoods import Gaussian
|
||||
from matplotlib import pyplot as plt
|
||||
|
||||
import numpy as np
|
||||
_np.random.seed(0)
|
||||
data = GPy.util.datasets.oil()
|
||||
try:
|
||||
import pods
|
||||
data = pods.datasets.oil()
|
||||
except ImportError:
|
||||
data = GPy.util.datasets.oil()
|
||||
|
||||
kernel = GPy.kern.rbf_inv(Q, 1., [.1] * Q, ARD=True) + GPy.kern.bias(Q, _np.exp(-2))
|
||||
|
||||
kernel = GPy.kern.RBF(Q, 1., 1. / _np.random.uniform(0, 1, (Q,)), ARD=True) # + GPy.kern.Bias(Q, _np.exp(-2))
|
||||
Y = data['X'][:N]
|
||||
Yn = Gaussian(Y, normalize=True)
|
||||
m = GPy.models.BayesianGPLVM(Yn, Q, kernel=kernel, num_inducing=num_inducing, **k)
|
||||
m = GPy.models.BayesianGPLVM(Y, Q, kernel=kernel, num_inducing=num_inducing, **k)
|
||||
m.data_labels = data['Y'][:N].argmax(axis=1)
|
||||
m['noise'] = Yn.Y.var() / 100.
|
||||
|
||||
if optimize:
|
||||
m.optimize('scg', messages=verbose, max_iters=max_iters, gtol=.05)
|
||||
m.optimize('bfgs', messages=verbose, max_iters=max_iters, gtol=.05)
|
||||
|
||||
if plot:
|
||||
y = m.likelihood.Y[0, :]
|
||||
fig, (latent_axes, sense_axes) = plt.subplots(1, 2)
|
||||
m.plot_latent(ax=latent_axes)
|
||||
data_show = GPy.util.visualize.vector_show(y)
|
||||
lvm_visualizer = GPy.util.visualize.lvm_dimselect(m.X[0, :], # @UnusedVariable
|
||||
m, data_show, latent_axes=latent_axes, sense_axes=sense_axes)
|
||||
m.plot_latent(ax=latent_axes, labels=m.data_labels)
|
||||
data_show = GPy.plotting.matplot_dep.visualize.vector_show((m.Y[0, :]))
|
||||
lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm_dimselect(m.X.mean.values[0:1, :], # @UnusedVariable
|
||||
m, data_show, latent_axes=latent_axes, sense_axes=sense_axes, labels=m.data_labels)
|
||||
raw_input('Press enter to finish')
|
||||
plt.close(fig)
|
||||
return m
|
||||
|
||||
def _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim=False):
|
||||
x = _np.linspace(0, 4 * _np.pi, N)[:, None]
|
||||
s1 = _np.vectorize(lambda x: _np.sin(x))
|
||||
s2 = _np.vectorize(lambda x: _np.cos(x))
|
||||
s3 = _np.vectorize(lambda x:-_np.exp(-_np.cos(2 * x)))
|
||||
sS = _np.vectorize(lambda x: _np.sin(2 * x))
|
||||
def ssgplvm_oil(optimize=True, verbose=1, plot=True, N=200, Q=7, num_inducing=40, max_iters=1000, **k):
|
||||
import GPy
|
||||
from matplotlib import pyplot as plt
|
||||
import pods
|
||||
|
||||
s1 = s1(x)
|
||||
s2 = s2(x)
|
||||
s3 = s3(x)
|
||||
sS = sS(x)
|
||||
_np.random.seed(0)
|
||||
data = pods.datasets.oil()
|
||||
|
||||
S1 = _np.hstack([s1, sS])
|
||||
S2 = _np.hstack([s2, s3, sS])
|
||||
S3 = _np.hstack([s3, sS])
|
||||
kernel = GPy.kern.RBF(Q, 1., 1. / _np.random.uniform(0, 1, (Q,)), ARD=True) # + GPy.kern.Bias(Q, _np.exp(-2))
|
||||
Y = data['X'][:N]
|
||||
m = GPy.models.SSGPLVM(Y, Q, kernel=kernel, num_inducing=num_inducing, **k)
|
||||
m.data_labels = data['Y'][:N].argmax(axis=1)
|
||||
|
||||
Y1 = S1.dot(_np.random.randn(S1.shape[1], D1))
|
||||
Y2 = S2.dot(_np.random.randn(S2.shape[1], D2))
|
||||
Y3 = S3.dot(_np.random.randn(S3.shape[1], D3))
|
||||
if optimize:
|
||||
m.optimize('bfgs', messages=verbose, max_iters=max_iters, gtol=.05)
|
||||
|
||||
Y1 += .3 * _np.random.randn(*Y1.shape)
|
||||
Y2 += .2 * _np.random.randn(*Y2.shape)
|
||||
Y3 += .25 * _np.random.randn(*Y3.shape)
|
||||
if plot:
|
||||
fig, (latent_axes, sense_axes) = plt.subplots(1, 2)
|
||||
m.plot_latent(ax=latent_axes, labels=m.data_labels)
|
||||
data_show = GPy.plotting.matplot_dep.visualize.vector_show((m.Y[0, :]))
|
||||
lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm_dimselect(m.X.mean.values[0:1, :], # @UnusedVariable
|
||||
m, data_show, latent_axes=latent_axes, sense_axes=sense_axes, labels=m.data_labels)
|
||||
raw_input('Press enter to finish')
|
||||
plt.close(fig)
|
||||
return m
|
||||
|
||||
Y1 -= Y1.mean(0)
|
||||
Y2 -= Y2.mean(0)
|
||||
Y3 -= Y3.mean(0)
|
||||
Y1 /= Y1.std(0)
|
||||
Y2 /= Y2.std(0)
|
||||
Y3 /= Y3.std(0)
|
||||
def _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim=False):
|
||||
Q_signal = 4
|
||||
import GPy
|
||||
import numpy as np
|
||||
np.random.seed(3000)
|
||||
|
||||
k = GPy.kern.Matern32(Q_signal, 1., lengthscale=(np.random.uniform(1, 6, Q_signal)), ARD=1)
|
||||
for i in range(Q_signal):
|
||||
k += GPy.kern.PeriodicExponential(1, variance=1., active_dims=[i], period=3., lower=-2, upper=6)
|
||||
t = np.c_[[np.linspace(-1, 5, N) for _ in range(Q_signal)]].T
|
||||
K = k.K(t)
|
||||
s2, s1, s3, sS = np.random.multivariate_normal(np.zeros(K.shape[0]), K, size=(4))[:, :, None]
|
||||
|
||||
Y1, Y2, Y3, S1, S2, S3 = _generate_high_dimensional_output(D1, D2, D3, s1, s2, s3, sS)
|
||||
|
||||
slist = [sS, s1, s2, s3]
|
||||
slist_names = ["sS", "s1", "s2", "s3"]
|
||||
Ylist = [Y1, Y2, Y3]
|
||||
|
||||
if plot_sim:
|
||||
import pylab
|
||||
from matplotlib import pyplot as plt
|
||||
import matplotlib.cm as cm
|
||||
import itertools
|
||||
fig = pylab.figure("MRD Simulation Data", figsize=(8, 6))
|
||||
fig = plt.figure("MRD Simulation Data", figsize=(8, 6))
|
||||
fig.clf()
|
||||
ax = fig.add_subplot(2, 1, 1)
|
||||
labls = slist_names
|
||||
|
|
@ -229,30 +246,75 @@ def _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim=False):
|
|||
ax.legend()
|
||||
for i, Y in enumerate(Ylist):
|
||||
ax = fig.add_subplot(2, len(Ylist), len(Ylist) + 1 + i)
|
||||
ax.imshow(Y, aspect='auto', cmap=cm.gray) # @UndefinedVariable
|
||||
ax.imshow(Y, aspect='auto', cmap=cm.gray) # @UndefinedVariable
|
||||
ax.set_title("Y{}".format(i + 1))
|
||||
pylab.draw()
|
||||
pylab.tight_layout()
|
||||
plt.draw()
|
||||
plt.tight_layout()
|
||||
|
||||
return slist, [S1, S2, S3], Ylist
|
||||
|
||||
# def bgplvm_simulation_matlab_compare():
|
||||
# from GPy.util.datasets import simulation_BGPLVM
|
||||
# from GPy import kern
|
||||
# from GPy.models import BayesianGPLVM
|
||||
#
|
||||
# sim_data = simulation_BGPLVM()
|
||||
# Y = sim_data['Y']
|
||||
# mu = sim_data['mu']
|
||||
# num_inducing, [_, Q] = 3, mu.shape
|
||||
#
|
||||
# k = kern.linear(Q, ARD=True) + kern.bias(Q, _np.exp(-2)) + kern.white(Q, _np.exp(-2))
|
||||
# m = BayesianGPLVM(Y, Q, init="PCA", num_inducing=num_inducing, kernel=k,
|
||||
# _debug=False)
|
||||
# m.auto_scale_factor = True
|
||||
# m['noise'] = Y.var() / 100.
|
||||
# m['linear_variance'] = .01
|
||||
# return m
|
||||
def _simulate_sincos(D1, D2, D3, N, num_inducing, plot_sim=False):
|
||||
_np.random.seed(1234)
|
||||
|
||||
x = _np.linspace(0, 4 * _np.pi, N)[:, None]
|
||||
s1 = _np.vectorize(lambda x: _np.sin(x))
|
||||
s2 = _np.vectorize(lambda x: _np.cos(x) ** 2)
|
||||
s3 = _np.vectorize(lambda x:-_np.exp(-_np.cos(2 * x)))
|
||||
sS = _np.vectorize(lambda x: _np.cos(x))
|
||||
|
||||
s1 = s1(x)
|
||||
s2 = s2(x)
|
||||
s3 = s3(x)
|
||||
sS = sS(x)
|
||||
|
||||
s1 -= s1.mean(); s1 /= s1.std(0)
|
||||
s2 -= s2.mean(); s2 /= s2.std(0)
|
||||
s3 -= s3.mean(); s3 /= s3.std(0)
|
||||
sS -= sS.mean(); sS /= sS.std(0)
|
||||
|
||||
Y1, Y2, Y3, S1, S2, S3 = _generate_high_dimensional_output(D1, D2, D3, s1, s2, s3, sS)
|
||||
|
||||
slist = [sS, s1, s2, s3]
|
||||
slist_names = ["sS", "s1", "s2", "s3"]
|
||||
Ylist = [Y1, Y2, Y3]
|
||||
|
||||
if plot_sim:
|
||||
from matplotlib import pyplot as plt
|
||||
import matplotlib.cm as cm
|
||||
import itertools
|
||||
fig = plt.figure("MRD Simulation Data", figsize=(8, 6))
|
||||
fig.clf()
|
||||
ax = fig.add_subplot(2, 1, 1)
|
||||
labls = slist_names
|
||||
for S, lab in itertools.izip(slist, labls):
|
||||
ax.plot(S, label=lab)
|
||||
ax.legend()
|
||||
for i, Y in enumerate(Ylist):
|
||||
ax = fig.add_subplot(2, len(Ylist), len(Ylist) + 1 + i)
|
||||
ax.imshow(Y, aspect='auto', cmap=cm.gray) # @UndefinedVariable
|
||||
ax.set_title("Y{}".format(i + 1))
|
||||
plt.draw()
|
||||
plt.tight_layout()
|
||||
|
||||
return slist, [S1, S2, S3], Ylist
|
||||
|
||||
def _generate_high_dimensional_output(D1, D2, D3, s1, s2, s3, sS):
|
||||
S1 = _np.hstack([s1, sS])
|
||||
S2 = _np.hstack([s2, s3, sS])
|
||||
S3 = _np.hstack([s3, sS])
|
||||
Y1 = S1.dot(_np.random.randn(S1.shape[1], D1))
|
||||
Y2 = S2.dot(_np.random.randn(S2.shape[1], D2))
|
||||
Y3 = S3.dot(_np.random.randn(S3.shape[1], D3))
|
||||
Y1 += .3 * _np.random.randn(*Y1.shape)
|
||||
Y2 += .2 * _np.random.randn(*Y2.shape)
|
||||
Y3 += .25 * _np.random.randn(*Y3.shape)
|
||||
Y1 -= Y1.mean(0)
|
||||
Y2 -= Y2.mean(0)
|
||||
Y3 -= Y3.mean(0)
|
||||
Y1 /= Y1.std(0)
|
||||
Y2 /= Y2.std(0)
|
||||
Y3 /= Y3.std(0)
|
||||
return Y1, Y2, Y3, S1, S2, S3
|
||||
|
||||
def bgplvm_simulation(optimize=True, verbose=1,
|
||||
plot=True, plot_sim=False,
|
||||
|
|
@ -261,95 +323,181 @@ def bgplvm_simulation(optimize=True, verbose=1,
|
|||
from GPy import kern
|
||||
from GPy.models import BayesianGPLVM
|
||||
|
||||
D1, D2, D3, N, num_inducing, Q = 49, 30, 10, 12, 3, 10
|
||||
_, _, Ylist = _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim)
|
||||
D1, D2, D3, N, num_inducing, Q = 13, 5, 8, 45, 3, 9
|
||||
_, _, Ylist = _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim)
|
||||
Y = Ylist[0]
|
||||
k = kern.linear(Q, ARD=True)
|
||||
k = kern.Linear(Q, ARD=True) # + kern.white(Q, _np.exp(-2)) # + kern.bias(Q)
|
||||
# k = kern.RBF(Q, ARD=True, lengthscale=10.)
|
||||
m = BayesianGPLVM(Y, Q, init="PCA", num_inducing=num_inducing, kernel=k)
|
||||
m.X_variance = m.X_variance * .7
|
||||
m['noise'] = Y.var() / 100.
|
||||
m.X.variance[:] = _np.random.uniform(0, .01, m.X.shape)
|
||||
m.likelihood.variance = .1
|
||||
|
||||
if optimize:
|
||||
print "Optimizing model:"
|
||||
m.optimize('bfgs', messages=verbose, max_iters=max_iters,
|
||||
gtol=.05)
|
||||
if plot:
|
||||
m.X.plot("BGPLVM Latent Space 1D")
|
||||
m.kern.plot_ARD('BGPLVM Simulation ARD Parameters')
|
||||
return m
|
||||
|
||||
def ssgplvm_simulation(optimize=True, verbose=1,
|
||||
plot=True, plot_sim=False,
|
||||
max_iters=2e4, useGPU=False
|
||||
):
|
||||
from GPy import kern
|
||||
from GPy.models import SSGPLVM
|
||||
|
||||
D1, D2, D3, N, num_inducing, Q = 13, 5, 8, 45, 3, 9
|
||||
_, _, Ylist = _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim)
|
||||
Y = Ylist[0]
|
||||
k = kern.Linear(Q, ARD=True) # + kern.white(Q, _np.exp(-2)) # + kern.bias(Q)
|
||||
# k = kern.RBF(Q, ARD=True, lengthscale=10.)
|
||||
m = SSGPLVM(Y, Q, init="pca", num_inducing=num_inducing, kernel=k)
|
||||
m.X.variance[:] = _np.random.uniform(0, .01, m.X.shape)
|
||||
m.likelihood.variance = .1
|
||||
|
||||
if optimize:
|
||||
print "Optimizing model:"
|
||||
m.optimize('scg', messages=verbose, max_iters=max_iters,
|
||||
gtol=.05)
|
||||
if plot:
|
||||
m.plot_X_1d("BGPLVM Latent Space 1D")
|
||||
m.X.plot("SSGPLVM Latent Space 1D")
|
||||
m.kern.plot_ARD('SSGPLVM Simulation ARD Parameters')
|
||||
return m
|
||||
|
||||
def bgplvm_simulation_missing_data(optimize=True, verbose=1,
|
||||
plot=True, plot_sim=False,
|
||||
max_iters=2e4, percent_missing=.1,
|
||||
):
|
||||
from GPy import kern
|
||||
from GPy.models.bayesian_gplvm_minibatch import BayesianGPLVMMiniBatch
|
||||
|
||||
D1, D2, D3, N, num_inducing, Q = 13, 5, 8, 400, 3, 4
|
||||
_, _, Ylist = _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim)
|
||||
Y = Ylist[0]
|
||||
k = kern.Linear(Q, ARD=True) # + kern.white(Q, _np.exp(-2)) # + kern.bias(Q)
|
||||
|
||||
inan = _np.random.binomial(1, percent_missing, size=Y.shape).astype(bool) # 80% missing data
|
||||
Ymissing = Y.copy()
|
||||
Ymissing[inan] = _np.nan
|
||||
|
||||
m = BayesianGPLVMMiniBatch(Ymissing, Q, init="random", num_inducing=num_inducing,
|
||||
kernel=k, missing_data=True)
|
||||
|
||||
m.Yreal = Y
|
||||
|
||||
if optimize:
|
||||
print "Optimizing model:"
|
||||
m.optimize('bfgs', messages=verbose, max_iters=max_iters,
|
||||
gtol=.05)
|
||||
if plot:
|
||||
m.X.plot("BGPLVM Latent Space 1D")
|
||||
m.kern.plot_ARD('BGPLVM Simulation ARD Parameters')
|
||||
return m
|
||||
|
||||
|
||||
def mrd_simulation(optimize=True, verbose=True, plot=True, plot_sim=True, **kw):
|
||||
from GPy import kern
|
||||
from GPy.models import MRD
|
||||
from GPy.likelihoods import Gaussian
|
||||
|
||||
D1, D2, D3, N, num_inducing, Q = 60, 20, 36, 60, 6, 5
|
||||
_, _, Ylist = _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim)
|
||||
likelihood_list = [Gaussian(x, normalize=True) for x in Ylist]
|
||||
_, _, Ylist = _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim)
|
||||
|
||||
k = kern.linear(Q, ARD=True)# + kern.bias(Q, _np.exp(-2)) + kern.white(Q, _np.exp(-2))
|
||||
m = MRD(likelihood_list, input_dim=Q, num_inducing=num_inducing, kernels=k, initx="", initz='permute', **kw)
|
||||
m.ensure_default_constraints()
|
||||
# Ylist = [Ylist[0]]
|
||||
k = kern.Linear(Q, ARD=True)
|
||||
m = MRD(Ylist, input_dim=Q, num_inducing=num_inducing, kernel=k, initx="PCA_concat", initz='permute', **kw)
|
||||
|
||||
m['.*noise'] = [Y.var() / 40. for Y in Ylist]
|
||||
|
||||
for i, bgplvm in enumerate(m.bgplvms):
|
||||
m['{}_noise'.format(i)] = 1 #bgplvm.likelihood.Y.var() / 500.
|
||||
bgplvm.X_variance = bgplvm.X_variance #* .1
|
||||
if optimize:
|
||||
print "Optimizing Model:"
|
||||
m.optimize(messages=verbose, max_iters=8e3, gtol=.1)
|
||||
m.optimize(messages=verbose, max_iters=8e3)
|
||||
if plot:
|
||||
m.plot_X_1d("MRD Latent Space 1D")
|
||||
m.X.plot("MRD Latent Space 1D")
|
||||
m.plot_scales("MRD Scales")
|
||||
return m
|
||||
|
||||
def mrd_simulation_missing_data(optimize=True, verbose=True, plot=True, plot_sim=True, **kw):
|
||||
from GPy import kern
|
||||
from GPy.models import MRD
|
||||
|
||||
D1, D2, D3, N, num_inducing, Q = 60, 20, 36, 60, 6, 5
|
||||
_, _, Ylist = _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim)
|
||||
|
||||
# Ylist = [Ylist[0]]
|
||||
k = kern.Linear(Q, ARD=True)
|
||||
inanlist = []
|
||||
|
||||
for Y in Ylist:
|
||||
inan = _np.random.binomial(1, .6, size=Y.shape).astype(bool)
|
||||
inanlist.append(inan)
|
||||
Y[inan] = _np.nan
|
||||
|
||||
m = MRD(Ylist, input_dim=Q, num_inducing=num_inducing,
|
||||
kernel=k, inference_method=None,
|
||||
initx="random", initz='permute', **kw)
|
||||
|
||||
if optimize:
|
||||
print "Optimizing Model:"
|
||||
m.optimize('bfgs', messages=verbose, max_iters=8e3, gtol=.1)
|
||||
if plot:
|
||||
m.X.plot("MRD Latent Space 1D")
|
||||
m.plot_scales("MRD Scales")
|
||||
return m
|
||||
|
||||
def brendan_faces(optimize=True, verbose=True, plot=True):
|
||||
import GPy
|
||||
import pods
|
||||
|
||||
data = GPy.util.datasets.brendan_faces()
|
||||
data = pods.datasets.brendan_faces()
|
||||
Q = 2
|
||||
Y = data['Y']
|
||||
Yn = Y - Y.mean()
|
||||
Yn /= Yn.std()
|
||||
|
||||
m = GPy.models.GPLVM(Yn, Q)
|
||||
m = GPy.models.BayesianGPLVM(Yn, Q, num_inducing=20)
|
||||
|
||||
# optimize
|
||||
m.constrain('rbf|noise|white', GPy.core.transformations.logexp_clipped())
|
||||
|
||||
if optimize: m.optimize('scg', messages=verbose, max_iters=1000)
|
||||
if optimize: m.optimize('bfgs', messages=verbose, max_iters=1000)
|
||||
|
||||
if plot:
|
||||
ax = m.plot_latent(which_indices=(0, 1))
|
||||
y = m.likelihood.Y[0, :]
|
||||
data_show = GPy.util.visualize.image_show(y[None, :], dimensions=(20, 28), transpose=True, order='F', invert=False, scale=False)
|
||||
GPy.util.visualize.lvm(m.X[0, :].copy(), m, data_show, ax)
|
||||
y = m.Y[0, :]
|
||||
data_show = GPy.plotting.matplot_dep.visualize.image_show(y[None, :], dimensions=(20, 28), transpose=True, order='F', invert=False, scale=False)
|
||||
lvm = GPy.plotting.matplot_dep.visualize.lvm(m.X.mean[0, :].copy(), m, data_show, ax)
|
||||
raw_input('Press enter to finish')
|
||||
|
||||
return m
|
||||
|
||||
def olivetti_faces(optimize=True, verbose=True, plot=True):
|
||||
import GPy
|
||||
import pods
|
||||
|
||||
data = GPy.util.datasets.olivetti_faces()
|
||||
data = pods.datasets.olivetti_faces()
|
||||
Q = 2
|
||||
Y = data['Y']
|
||||
Yn = Y - Y.mean()
|
||||
Yn /= Yn.std()
|
||||
|
||||
m = GPy.models.GPLVM(Yn, Q)
|
||||
if optimize: m.optimize('scg', messages=verbose, max_iters=1000)
|
||||
m = GPy.models.BayesianGPLVM(Yn, Q, num_inducing=20)
|
||||
|
||||
if optimize: m.optimize('bfgs', messages=verbose, max_iters=1000)
|
||||
if plot:
|
||||
ax = m.plot_latent(which_indices=(0, 1))
|
||||
y = m.likelihood.Y[0, :]
|
||||
data_show = GPy.util.visualize.image_show(y[None, :], dimensions=(112, 92), transpose=False, invert=False, scale=False)
|
||||
GPy.util.visualize.lvm(m.X[0, :].copy(), m, data_show, ax)
|
||||
y = m.Y[0, :]
|
||||
data_show = GPy.plotting.matplot_dep.visualize.image_show(y[None, :], dimensions=(112, 92), transpose=False, invert=False, scale=False)
|
||||
lvm = GPy.plotting.matplot_dep.visualize.lvm(m.X.mean[0, :].copy(), m, data_show, ax)
|
||||
raw_input('Press enter to finish')
|
||||
|
||||
return m
|
||||
|
||||
def stick_play(range=None, frame_rate=15, optimize=False, verbose=True, plot=True):
|
||||
import GPy
|
||||
data = GPy.util.datasets.osu_run1()
|
||||
import pods
|
||||
|
||||
data = pods.datasets.osu_run1()
|
||||
# optimize
|
||||
if range == None:
|
||||
Y = data['Y'].copy()
|
||||
|
|
@ -357,43 +505,46 @@ def stick_play(range=None, frame_rate=15, optimize=False, verbose=True, plot=Tru
|
|||
Y = data['Y'][range[0]:range[1], :].copy()
|
||||
if plot:
|
||||
y = Y[0, :]
|
||||
data_show = GPy.util.visualize.stick_show(y[None, :], connect=data['connect'])
|
||||
GPy.util.visualize.data_play(Y, data_show, frame_rate)
|
||||
data_show = GPy.plotting.matplot_dep.visualize.stick_show(y[None, :], connect=data['connect'])
|
||||
GPy.plotting.matplot_dep.visualize.data_play(Y, data_show, frame_rate)
|
||||
return Y
|
||||
|
||||
def stick(kernel=None, optimize=True, verbose=True, plot=True):
|
||||
from matplotlib import pyplot as plt
|
||||
import GPy
|
||||
import pods
|
||||
|
||||
data = GPy.util.datasets.osu_run1()
|
||||
data = pods.datasets.osu_run1()
|
||||
# optimize
|
||||
m = GPy.models.GPLVM(data['Y'], 2, kernel=kernel)
|
||||
if optimize: m.optimize(messages=verbose, max_f_eval=10000)
|
||||
if plot and GPy.util.visualize.visual_available:
|
||||
if optimize: m.optimize('bfgs', messages=verbose, max_f_eval=10000)
|
||||
if plot:
|
||||
plt.clf
|
||||
ax = m.plot_latent()
|
||||
y = m.likelihood.Y[0, :]
|
||||
data_show = GPy.util.visualize.stick_show(y[None, :], connect=data['connect'])
|
||||
GPy.util.visualize.lvm(m.X[0, :].copy(), m, data_show, ax)
|
||||
y = m.Y[0, :]
|
||||
data_show = GPy.plotting.matplot_dep.visualize.stick_show(y[None, :], connect=data['connect'])
|
||||
lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm(m.X[:1, :].copy(), m, data_show, latent_axes=ax)
|
||||
raw_input('Press enter to finish')
|
||||
|
||||
lvm_visualizer.close()
|
||||
data_show.close()
|
||||
return m
|
||||
|
||||
def bcgplvm_linear_stick(kernel=None, optimize=True, verbose=True, plot=True):
|
||||
from matplotlib import pyplot as plt
|
||||
import GPy
|
||||
import pods
|
||||
|
||||
data = GPy.util.datasets.osu_run1()
|
||||
data = pods.datasets.osu_run1()
|
||||
# optimize
|
||||
mapping = GPy.mappings.Linear(data['Y'].shape[1], 2)
|
||||
m = GPy.models.BCGPLVM(data['Y'], 2, kernel=kernel, mapping=mapping)
|
||||
if optimize: m.optimize(messages=verbose, max_f_eval=10000)
|
||||
if plot and GPy.util.visualize.visual_available:
|
||||
if plot and GPy.plotting.matplot_dep.visualize.visual_available:
|
||||
plt.clf
|
||||
ax = m.plot_latent()
|
||||
y = m.likelihood.Y[0, :]
|
||||
data_show = GPy.util.visualize.stick_show(y[None, :], connect=data['connect'])
|
||||
GPy.util.visualize.lvm(m.X[0, :].copy(), m, data_show, ax)
|
||||
data_show = GPy.plotting.matplot_dep.visualize.stick_show(y[None, :], connect=data['connect'])
|
||||
GPy.plotting.matplot_dep.visualize.lvm(m.X[0, :].copy(), m, data_show, ax)
|
||||
raw_input('Press enter to finish')
|
||||
|
||||
return m
|
||||
|
|
@ -401,32 +552,33 @@ def bcgplvm_linear_stick(kernel=None, optimize=True, verbose=True, plot=True):
|
|||
def bcgplvm_stick(kernel=None, optimize=True, verbose=True, plot=True):
|
||||
from matplotlib import pyplot as plt
|
||||
import GPy
|
||||
import pods
|
||||
|
||||
data = GPy.util.datasets.osu_run1()
|
||||
data = pods.datasets.osu_run1()
|
||||
# optimize
|
||||
back_kernel=GPy.kern.rbf(data['Y'].shape[1], lengthscale=5.)
|
||||
back_kernel = GPy.kern.RBF(data['Y'].shape[1], lengthscale=5.)
|
||||
mapping = GPy.mappings.Kernel(X=data['Y'], output_dim=2, kernel=back_kernel)
|
||||
m = GPy.models.BCGPLVM(data['Y'], 2, kernel=kernel, mapping=mapping)
|
||||
if optimize: m.optimize(messages=verbose, max_f_eval=10000)
|
||||
if plot and GPy.util.visualize.visual_available:
|
||||
if plot and GPy.plotting.matplot_dep.visualize.visual_available:
|
||||
plt.clf
|
||||
ax = m.plot_latent()
|
||||
y = m.likelihood.Y[0, :]
|
||||
data_show = GPy.util.visualize.stick_show(y[None, :], connect=data['connect'])
|
||||
GPy.util.visualize.lvm(m.X[0, :].copy(), m, data_show, ax)
|
||||
raw_input('Press enter to finish')
|
||||
data_show = GPy.plotting.matplot_dep.visualize.stick_show(y[None, :], connect=data['connect'])
|
||||
GPy.plotting.matplot_dep.visualize.lvm(m.X[0, :].copy(), m, data_show, ax)
|
||||
# raw_input('Press enter to finish')
|
||||
|
||||
return m
|
||||
|
||||
def robot_wireless(optimize=True, verbose=True, plot=True):
|
||||
from matplotlib import pyplot as plt
|
||||
import GPy
|
||||
import pods
|
||||
|
||||
data = GPy.util.datasets.robot_wireless()
|
||||
data = pods.datasets.robot_wireless()
|
||||
# optimize
|
||||
m = GPy.models.GPLVM(data['Y'], 2)
|
||||
m = GPy.models.BayesianGPLVM(data['Y'], 4, num_inducing=25)
|
||||
if optimize: m.optimize(messages=verbose, max_f_eval=10000)
|
||||
m._set_params(m._get_params())
|
||||
if plot:
|
||||
m.plot_latent()
|
||||
|
||||
|
|
@ -435,23 +587,33 @@ def robot_wireless(optimize=True, verbose=True, plot=True):
|
|||
def stick_bgplvm(model=None, optimize=True, verbose=True, plot=True):
|
||||
from GPy.models import BayesianGPLVM
|
||||
from matplotlib import pyplot as plt
|
||||
import numpy as np
|
||||
import GPy
|
||||
import pods
|
||||
|
||||
data = GPy.util.datasets.osu_run1()
|
||||
data = pods.datasets.osu_run1()
|
||||
Q = 6
|
||||
kernel = GPy.kern.rbf(Q, ARD=True) + GPy.kern.bias(Q, _np.exp(-2)) + GPy.kern.white(Q, _np.exp(-2))
|
||||
kernel = GPy.kern.RBF(Q, lengthscale=np.repeat(.5, Q), ARD=True)
|
||||
m = BayesianGPLVM(data['Y'], Q, init="PCA", num_inducing=20, kernel=kernel)
|
||||
|
||||
m.data = data
|
||||
m.likelihood.variance = 0.001
|
||||
|
||||
# optimize
|
||||
m.ensure_default_constraints()
|
||||
if optimize: m.optimize('scg', messages=verbose, max_iters=200, xtol=1e-300, ftol=1e-300)
|
||||
m._set_params(m._get_params())
|
||||
try:
|
||||
if optimize: m.optimize('bfgs', messages=verbose, max_iters=5e3, bfgs_factor=10)
|
||||
except KeyboardInterrupt:
|
||||
print "Keyboard interrupt, continuing to plot and return"
|
||||
|
||||
if plot:
|
||||
plt.clf, (latent_axes, sense_axes) = plt.subplots(1, 2)
|
||||
fig, (latent_axes, sense_axes) = plt.subplots(1, 2)
|
||||
plt.sca(latent_axes)
|
||||
m.plot_latent()
|
||||
y = m.likelihood.Y[0, :].copy()
|
||||
data_show = GPy.util.visualize.stick_show(y[None, :], connect=data['connect'])
|
||||
GPy.util.visualize.lvm_dimselect(m.X[0, :].copy(), m, data_show, latent_axes=latent_axes, sense_axes=sense_axes)
|
||||
m.plot_latent(ax=latent_axes)
|
||||
y = m.Y[:1, :].copy()
|
||||
data_show = GPy.plotting.matplot_dep.visualize.stick_show(y, connect=data['connect'])
|
||||
dim_select = GPy.plotting.matplot_dep.visualize.lvm_dimselect(m.X.mean[:1, :].copy(), m, data_show, latent_axes=latent_axes, sense_axes=sense_axes)
|
||||
fig.canvas.draw()
|
||||
fig.canvas.show()
|
||||
raw_input('Press enter to finish')
|
||||
|
||||
return m
|
||||
|
|
@ -459,20 +621,50 @@ def stick_bgplvm(model=None, optimize=True, verbose=True, plot=True):
|
|||
|
||||
def cmu_mocap(subject='35', motion=['01'], in_place=True, optimize=True, verbose=True, plot=True):
|
||||
import GPy
|
||||
import pods
|
||||
|
||||
data = GPy.util.datasets.cmu_mocap(subject, motion)
|
||||
data = pods.datasets.cmu_mocap(subject, motion)
|
||||
if in_place:
|
||||
# Make figure move in place.
|
||||
data['Y'][:, 0:3] = 0.0
|
||||
m = GPy.models.GPLVM(data['Y'], 2, normalize_Y=True)
|
||||
Y = data['Y']
|
||||
Y_mean = Y.mean(0)
|
||||
Y_std = Y.std(0)
|
||||
m = GPy.models.GPLVM((Y - Y_mean) / Y_std, 2)
|
||||
|
||||
if optimize: m.optimize(messages=verbose, max_f_eval=10000)
|
||||
if plot:
|
||||
ax = m.plot_latent()
|
||||
y = m.likelihood.Y[0, :]
|
||||
data_show = GPy.util.visualize.skeleton_show(y[None, :], data['skel'])
|
||||
lvm_visualizer = GPy.util.visualize.lvm(m.X[0, :].copy(), m, data_show, ax)
|
||||
y = m.Y[0, :]
|
||||
data_show = GPy.plotting.matplot_dep.visualize.skeleton_show(y[None, :], data['skel'])
|
||||
lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm(m.X[0].copy(), m, data_show, latent_axes=ax)
|
||||
raw_input('Press enter to finish')
|
||||
lvm_visualizer.close()
|
||||
data_show.close()
|
||||
|
||||
return m
|
||||
|
||||
def ssgplvm_simulation_linear():
|
||||
import numpy as np
|
||||
import GPy
|
||||
N, D, Q = 1000, 20, 5
|
||||
pi = 0.2
|
||||
|
||||
def sample_X(Q, pi):
|
||||
x = np.empty(Q)
|
||||
dies = np.random.rand(Q)
|
||||
for q in xrange(Q):
|
||||
if dies[q] < pi:
|
||||
x[q] = np.random.randn()
|
||||
else:
|
||||
x[q] = 0.
|
||||
return x
|
||||
|
||||
Y = np.empty((N, D))
|
||||
X = np.empty((N, Q))
|
||||
# Generate data from random sampled weight matrices
|
||||
for n in xrange(N):
|
||||
X[n] = sample_X(Q, pi)
|
||||
w = np.random.randn(D, Q)
|
||||
Y[n] = np.dot(w, X[n])
|
||||
|
||||
|
|
|
|||
|
|
@ -1,7 +1,13 @@
|
|||
# Copyright (c) 2014, Alan Saul
|
||||
# Licensed under the BSD 3-clause license (see LICENSE.txt)
|
||||
|
||||
import GPy
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from GPy.util import datasets
|
||||
try:
|
||||
import matplotlib.pyplot as plt
|
||||
except:
|
||||
pass
|
||||
|
||||
def student_t_approx(optimize=True, plot=True):
|
||||
"""
|
||||
|
|
@ -30,47 +36,53 @@ def student_t_approx(optimize=True, plot=True):
|
|||
#Yc = Yc/Yc.max()
|
||||
|
||||
#Add student t random noise to datapoints
|
||||
deg_free = 5
|
||||
deg_free = 1
|
||||
print "Real noise: ", real_std
|
||||
initial_var_guess = 0.5
|
||||
edited_real_sd = initial_var_guess
|
||||
|
||||
# Kernel object
|
||||
kernel1 = GPy.kern.rbf(X.shape[1]) + GPy.kern.white(X.shape[1])
|
||||
kernel2 = kernel1.copy()
|
||||
kernel3 = kernel1.copy()
|
||||
kernel4 = kernel1.copy()
|
||||
kernel1 = GPy.kern.RBF(X.shape[1]) + GPy.kern.White(X.shape[1])
|
||||
kernel2 = GPy.kern.RBF(X.shape[1]) + GPy.kern.White(X.shape[1])
|
||||
kernel3 = GPy.kern.RBF(X.shape[1]) + GPy.kern.White(X.shape[1])
|
||||
kernel4 = GPy.kern.RBF(X.shape[1]) + GPy.kern.White(X.shape[1])
|
||||
|
||||
#Gaussian GP model on clean data
|
||||
m1 = GPy.models.GPRegression(X, Y.copy(), kernel=kernel1)
|
||||
# optimize
|
||||
m1.ensure_default_constraints()
|
||||
m1.constrain_fixed('white', 1e-5)
|
||||
m1['.*white'].constrain_fixed(1e-5)
|
||||
m1.randomize()
|
||||
|
||||
#Gaussian GP model on corrupt data
|
||||
m2 = GPy.models.GPRegression(X, Yc.copy(), kernel=kernel2)
|
||||
m2.ensure_default_constraints()
|
||||
m2.constrain_fixed('white', 1e-5)
|
||||
m2['.*white'].constrain_fixed(1e-5)
|
||||
m2.randomize()
|
||||
|
||||
#Student t GP model on clean data
|
||||
t_distribution = GPy.likelihoods.noise_model_constructors.student_t(deg_free=deg_free, sigma2=edited_real_sd)
|
||||
stu_t_likelihood = GPy.likelihoods.Laplace(Y.copy(), t_distribution)
|
||||
m3 = GPy.models.GPRegression(X, Y.copy(), kernel3, likelihood=stu_t_likelihood)
|
||||
m3.ensure_default_constraints()
|
||||
m3.constrain_bounded('t_noise', 1e-6, 10.)
|
||||
m3.constrain_fixed('white', 1e-5)
|
||||
t_distribution = GPy.likelihoods.StudentT(deg_free=deg_free, sigma2=edited_real_sd)
|
||||
laplace_inf = GPy.inference.latent_function_inference.Laplace()
|
||||
m3 = GPy.core.GP(X, Y.copy(), kernel3, likelihood=t_distribution, inference_method=laplace_inf)
|
||||
m3['.*t_scale2'].constrain_bounded(1e-6, 10.)
|
||||
m3['.*white'].constrain_fixed(1e-5)
|
||||
m3.randomize()
|
||||
|
||||
#Student t GP model on corrupt data
|
||||
t_distribution = GPy.likelihoods.noise_model_constructors.student_t(deg_free=deg_free, sigma2=edited_real_sd)
|
||||
corrupt_stu_t_likelihood = GPy.likelihoods.Laplace(Yc.copy(), t_distribution)
|
||||
m4 = GPy.models.GPRegression(X, Yc.copy(), kernel4, likelihood=corrupt_stu_t_likelihood)
|
||||
m4.ensure_default_constraints()
|
||||
m4.constrain_bounded('t_noise', 1e-6, 10.)
|
||||
m4.constrain_fixed('white', 1e-5)
|
||||
t_distribution = GPy.likelihoods.StudentT(deg_free=deg_free, sigma2=edited_real_sd)
|
||||
laplace_inf = GPy.inference.latent_function_inference.Laplace()
|
||||
m4 = GPy.core.GP(X, Yc.copy(), kernel4, likelihood=t_distribution, inference_method=laplace_inf)
|
||||
m4['.*t_scale2'].constrain_bounded(1e-6, 10.)
|
||||
m4['.*white'].constrain_fixed(1e-5)
|
||||
m4.randomize()
|
||||
print m4
|
||||
debug=True
|
||||
if debug:
|
||||
m4.optimize(messages=1)
|
||||
import pylab as pb
|
||||
pb.plot(m4.X, m4.inference_method.f_hat)
|
||||
pb.plot(m4.X, m4.Y, 'rx')
|
||||
m4.plot()
|
||||
print m4
|
||||
return m4
|
||||
|
||||
if optimize:
|
||||
optimizer='scg'
|
||||
|
|
@ -115,6 +127,7 @@ def student_t_approx(optimize=True, plot=True):
|
|||
return m1, m2, m3, m4
|
||||
|
||||
def boston_example(optimize=True, plot=True):
|
||||
raise NotImplementedError("Needs updating")
|
||||
import sklearn
|
||||
from sklearn.cross_validation import KFold
|
||||
optimizer='bfgs'
|
||||
|
|
@ -143,8 +156,8 @@ def boston_example(optimize=True, plot=True):
|
|||
noise = 1e-1 #np.exp(-2)
|
||||
rbf_len = 0.5
|
||||
data_axis_plot = 4
|
||||
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])
|
||||
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])
|
||||
|
||||
#Baseline
|
||||
score_folds[0, n] = rmse(Y_test, np.mean(Y_train))
|
||||
|
|
@ -152,10 +165,9 @@ def boston_example(optimize=True, plot=True):
|
|||
#Gaussian GP
|
||||
print "Gauss GP"
|
||||
mgp = GPy.models.GPRegression(X_train.copy(), Y_train.copy(), kernel=kernelgp.copy())
|
||||
mgp.ensure_default_constraints()
|
||||
mgp.constrain_fixed('white', 1e-5)
|
||||
mgp['rbf_len'] = rbf_len
|
||||
mgp['noise'] = noise
|
||||
mgp.constrain_fixed('.*white', 1e-5)
|
||||
mgp['.*len'] = rbf_len
|
||||
mgp['.*noise'] = noise
|
||||
print mgp
|
||||
if optimize:
|
||||
mgp.optimize(optimizer=optimizer, messages=messages)
|
||||
|
|
@ -170,9 +182,8 @@ def boston_example(optimize=True, plot=True):
|
|||
g_distribution = GPy.likelihoods.noise_model_constructors.gaussian(variance=noise, N=N, D=D)
|
||||
g_likelihood = GPy.likelihoods.Laplace(Y_train.copy(), g_distribution)
|
||||
mg = GPy.models.GPRegression(X_train.copy(), Y_train.copy(), kernel=kernelstu.copy(), likelihood=g_likelihood)
|
||||
mg.ensure_default_constraints()
|
||||
mg.constrain_positive('noise_variance')
|
||||
mg.constrain_fixed('white', 1e-5)
|
||||
mg.constrain_fixed('.*white', 1e-5)
|
||||
mg['rbf_len'] = rbf_len
|
||||
mg['noise'] = noise
|
||||
print mg
|
||||
|
|
@ -190,11 +201,10 @@ def boston_example(optimize=True, plot=True):
|
|||
t_distribution = GPy.likelihoods.noise_model_constructors.student_t(deg_free=df, sigma2=noise)
|
||||
stu_t_likelihood = GPy.likelihoods.Laplace(Y_train.copy(), t_distribution)
|
||||
mstu_t = GPy.models.GPRegression(X_train.copy(), Y_train.copy(), kernel=kernelstu.copy(), likelihood=stu_t_likelihood)
|
||||
mstu_t.ensure_default_constraints()
|
||||
mstu_t.constrain_fixed('white', 1e-5)
|
||||
mstu_t.constrain_bounded('t_noise', 0.0001, 1000)
|
||||
mstu_t.constrain_fixed('.*white', 1e-5)
|
||||
mstu_t.constrain_bounded('.*t_scale2', 0.0001, 1000)
|
||||
mstu_t['rbf_len'] = rbf_len
|
||||
mstu_t['t_noise'] = noise
|
||||
mstu_t['.*t_scale2'] = noise
|
||||
print mstu_t
|
||||
if optimize:
|
||||
mstu_t.optimize(optimizer=optimizer, messages=messages)
|
||||
|
|
|
|||
|
|
@ -1,22 +1,29 @@
|
|||
# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
|
||||
# Copyright (c) 2012-2014, GPy authors (see AUTHORS.txt).
|
||||
# Licensed under the BSD 3-clause license (see LICENSE.txt)
|
||||
|
||||
"""
|
||||
Gaussian Processes regression examples
|
||||
"""
|
||||
import pylab as pb
|
||||
try:
|
||||
import pylab as pb
|
||||
except:
|
||||
pass
|
||||
import numpy as np
|
||||
import GPy
|
||||
|
||||
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()
|
||||
try:import pods
|
||||
except ImportError:
|
||||
print 'pods unavailable, see https://github.com/sods/ods for example datasets'
|
||||
return
|
||||
data = pods.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
|
||||
m.kern.lengthscale = 10.
|
||||
|
||||
if optimize:
|
||||
m.optimize('bfgs', max_iters=200)
|
||||
|
|
@ -25,79 +32,51 @@ def olympic_marathon_men(optimize=True, plot=True):
|
|||
|
||||
return m
|
||||
|
||||
def coregionalization_toy2(optimize=True, plot=True):
|
||||
def coregionalization_toy(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
|
||||
m = GPy.models.GPRegression(X, Y, kernel=k)
|
||||
m.constrain_fixed('.*rbf_var', 1.)
|
||||
m = GPy.models.GPCoregionalizedRegression(X_list=[X1,X2], Y_list=[Y1,Y2])
|
||||
|
||||
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())
|
||||
|
||||
slices = GPy.util.multioutput.get_slices([X1,X2])
|
||||
m.plot(fixed_inputs=[(1,0)],which_data_rows=slices[0],Y_metadata={'output_index':0})
|
||||
m.plot(fixed_inputs=[(1,1)],which_data_rows=slices[1],Y_metadata={'output_index':1},ax=pb.gca())
|
||||
return m
|
||||
|
||||
#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
|
||||
|
||||
def coregionalization_sparse(optimize=True, plot=True):
|
||||
"""
|
||||
A simple demonstration of coregionalization on two sinusoidal functions using sparse approximations.
|
||||
"""
|
||||
#fetch the data from the non sparse examples
|
||||
m = coregionalization_toy2(optimize=False, plot=False)
|
||||
X, Y = m.X, m.likelihood.Y
|
||||
#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
|
||||
|
||||
#construct a model
|
||||
m = GPy.models.SparseGPRegression(X,Y)
|
||||
m.constrain_fixed('iip_\d+_1') # don't optimize the inducing input indexes
|
||||
#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.
|
||||
|
||||
m = GPy.models.SparseGPCoregionalizedRegression(X_list=[X1,X2], Y_list=[Y1,Y2])
|
||||
|
||||
if optimize:
|
||||
m.optimize('bfgs', max_iters=100, messages=1)
|
||||
m.optimize('bfgs', max_iters=100)
|
||||
|
||||
if plot:
|
||||
m.plot(fixed_inputs=[(1,0)])
|
||||
m.plot(fixed_inputs=[(1,1)], ax=pb.gca())
|
||||
slices = GPy.util.multioutput.get_slices([X1,X2])
|
||||
m.plot(fixed_inputs=[(1,0)],which_data_rows=slices[0],Y_metadata={'output_index':0})
|
||||
m.plot(fixed_inputs=[(1,1)],which_data_rows=slices[1],Y_metadata={'output_index':1},ax=pb.gca())
|
||||
pb.ylim(-3,)
|
||||
|
||||
return m
|
||||
|
||||
|
|
@ -107,7 +86,11 @@ def epomeo_gpx(max_iters=200, optimize=True, plot=True):
|
|||
from the Mount Epomeo runs. Requires gpxpy to be installed on your system
|
||||
to load in the data.
|
||||
"""
|
||||
data = GPy.util.datasets.epomeo_gpx()
|
||||
try:import pods
|
||||
except ImportError:
|
||||
print 'pods unavailable, see https://github.com/sods/ods for example datasets'
|
||||
return
|
||||
data = pods.datasets.epomeo_gpx()
|
||||
num_data_list = []
|
||||
for Xpart in data['X']:
|
||||
num_data_list.append(Xpart.shape[0])
|
||||
|
|
@ -127,14 +110,14 @@ def epomeo_gpx(max_iters=200, optimize=True, plot=True):
|
|||
Z = np.hstack((np.linspace(t[:,0].min(), t[:, 0].max(), num_inducing)[:, None],
|
||||
np.random.randint(0, 4, num_inducing)[:, None]))
|
||||
|
||||
k1 = GPy.kern.rbf(1)
|
||||
k2 = GPy.kern.coregionalize(output_dim=5, rank=5)
|
||||
k1 = GPy.kern.RBF(1)
|
||||
k2 = GPy.kern.Coregionalize(output_dim=5, rank=5)
|
||||
k = k1**k2
|
||||
|
||||
m = GPy.models.SparseGPRegression(t, Y, kernel=k, Z=Z, normalize_Y=True)
|
||||
m.constrain_fixed('.*rbf_var', 1.)
|
||||
m.constrain_fixed('iip')
|
||||
m.constrain_bounded('noise_variance', 1e-3, 1e-1)
|
||||
m.constrain_fixed('.*variance', 1.)
|
||||
m.inducing_inputs.constrain_fixed()
|
||||
m.Gaussian_noise.variance.constrain_bounded(1e-3, 1e-1)
|
||||
m.optimize(max_iters=max_iters,messages=True)
|
||||
|
||||
return m
|
||||
|
|
@ -150,13 +133,17 @@ def multiple_optima(gene_number=937, resolution=80, model_restarts=10, seed=1000
|
|||
length_scales = np.linspace(0.1, 60., resolution)
|
||||
log_SNRs = np.linspace(-3., 4., resolution)
|
||||
|
||||
data = GPy.util.datasets.della_gatta_TRP63_gene_expression(data_set='della_gatta',gene_number=gene_number)
|
||||
try:import pods
|
||||
except ImportError:
|
||||
print 'pods unavailable, see https://github.com/sods/ods for example datasets'
|
||||
return
|
||||
data = pods.datasets.della_gatta_TRP63_gene_expression(data_set='della_gatta',gene_number=gene_number)
|
||||
# data['Y'] = data['Y'][0::2, :]
|
||||
# data['X'] = data['X'][0::2, :]
|
||||
|
||||
data['Y'] = data['Y'] - np.mean(data['Y'])
|
||||
|
||||
lls = GPy.examples.regression._contour_data(data, length_scales, log_SNRs, GPy.kern.rbf)
|
||||
lls = GPy.examples.regression._contour_data(data, length_scales, log_SNRs, GPy.kern.RBF)
|
||||
if plot:
|
||||
pb.contour(length_scales, log_SNRs, np.exp(lls), 20, cmap=pb.cm.jet)
|
||||
ax = pb.gca()
|
||||
|
|
@ -172,20 +159,20 @@ def multiple_optima(gene_number=937, resolution=80, model_restarts=10, seed=1000
|
|||
optim_point_y = np.empty(2)
|
||||
np.random.seed(seed=seed)
|
||||
for i in range(0, model_restarts):
|
||||
# kern = GPy.kern.rbf(1, variance=np.random.exponential(1.), lengthscale=np.random.exponential(50.))
|
||||
kern = GPy.kern.rbf(1, variance=np.random.uniform(1e-3, 1), lengthscale=np.random.uniform(5, 50))
|
||||
# kern = GPy.kern.RBF(1, variance=np.random.exponential(1.), lengthscale=np.random.exponential(50.))
|
||||
kern = GPy.kern.RBF(1, variance=np.random.uniform(1e-3, 1), lengthscale=np.random.uniform(5, 50))
|
||||
|
||||
m = GPy.models.GPRegression(data['X'], data['Y'], kernel=kern)
|
||||
m['noise_variance'] = np.random.uniform(1e-3, 1)
|
||||
optim_point_x[0] = m['rbf_lengthscale']
|
||||
optim_point_y[0] = np.log10(m['rbf_variance']) - np.log10(m['noise_variance']);
|
||||
m.likelihood.variance = np.random.uniform(1e-3, 1)
|
||||
optim_point_x[0] = m.rbf.lengthscale
|
||||
optim_point_y[0] = np.log10(m.rbf.variance) - np.log10(m.likelihood.variance);
|
||||
|
||||
# optimize
|
||||
if optimize:
|
||||
m.optimize('scg', xtol=1e-6, ftol=1e-6, max_iters=max_iters)
|
||||
|
||||
optim_point_x[1] = m['rbf_lengthscale']
|
||||
optim_point_y[1] = np.log10(m['rbf_variance']) - np.log10(m['noise_variance']);
|
||||
optim_point_x[1] = m.rbf.lengthscale
|
||||
optim_point_y[1] = np.log10(m.rbf.variance) - np.log10(m.likelihood.variance);
|
||||
|
||||
if plot:
|
||||
pb.arrow(optim_point_x[0], optim_point_y[0], optim_point_x[1] - optim_point_x[0], optim_point_y[1] - optim_point_y[0], label=str(i), head_length=1, head_width=0.5, fc='k', ec='k')
|
||||
|
|
@ -196,7 +183,7 @@ def multiple_optima(gene_number=937, resolution=80, model_restarts=10, seed=1000
|
|||
ax.set_ylim(ylim)
|
||||
return m # (models, lls)
|
||||
|
||||
def _contour_data(data, length_scales, log_SNRs, kernel_call=GPy.kern.rbf):
|
||||
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.
|
||||
|
|
@ -216,7 +203,7 @@ def _contour_data(data, length_scales, log_SNRs, kernel_call=GPy.kern.rbf):
|
|||
noise_var = total_var / (1. + SNR)
|
||||
signal_var = total_var - noise_var
|
||||
model.kern['.*variance'] = signal_var
|
||||
model['noise_variance'] = noise_var
|
||||
model.likelihood.variance = noise_var
|
||||
length_scale_lls = []
|
||||
|
||||
for length_scale in length_scales:
|
||||
|
|
@ -230,13 +217,17 @@ def _contour_data(data, length_scales, log_SNRs, kernel_call=GPy.kern.rbf):
|
|||
|
||||
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()
|
||||
try:import pods
|
||||
except ImportError:
|
||||
print 'pods unavailable, see https://github.com/sods/ods for example datasets'
|
||||
return
|
||||
data = pods.datasets.olympic_100m_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
|
||||
m.rbf.lengthscale = 10
|
||||
|
||||
if optimize:
|
||||
m.optimize('bfgs', max_iters=200)
|
||||
|
|
@ -247,7 +238,11 @@ def olympic_100m_men(optimize=True, plot=True):
|
|||
|
||||
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()
|
||||
try:import pods
|
||||
except ImportError:
|
||||
print 'pods unavailable, see https://github.com/sods/ods for example datasets'
|
||||
return
|
||||
data = pods.datasets.toy_rbf_1d()
|
||||
|
||||
# create simple GP Model
|
||||
m = GPy.models.GPRegression(data['X'], data['Y'])
|
||||
|
|
@ -261,7 +256,11 @@ def toy_rbf_1d(optimize=True, plot=True):
|
|||
|
||||
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()
|
||||
try:import pods
|
||||
except ImportError:
|
||||
print 'pods unavailable, see https://github.com/sods/ods for example datasets'
|
||||
return
|
||||
data = pods.datasets.toy_rbf_1d_50()
|
||||
|
||||
# create simple GP Model
|
||||
m = GPy.models.GPRegression(data['X'], data['Y'])
|
||||
|
|
@ -278,14 +277,15 @@ def toy_poisson_rbf_1d_laplace(optimize=True, plot=True):
|
|||
optimizer='scg'
|
||||
x_len = 30
|
||||
X = np.linspace(0, 10, x_len)[:, None]
|
||||
f_true = np.random.multivariate_normal(np.zeros(x_len), GPy.kern.rbf(1).K(X))
|
||||
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]
|
||||
|
||||
noise_model = GPy.likelihoods.poisson()
|
||||
likelihood = GPy.likelihoods.Laplace(Y,noise_model)
|
||||
kern = GPy.kern.RBF(1)
|
||||
poisson_lik = GPy.likelihoods.Poisson()
|
||||
laplace_inf = GPy.inference.latent_function_inference.Laplace()
|
||||
|
||||
# create simple GP Model
|
||||
m = GPy.models.GPRegression(X, Y, likelihood=likelihood)
|
||||
m = GPy.core.GP(X, Y, kernel=kern, likelihood=poisson_lik, inference_method=laplace_inf)
|
||||
|
||||
if optimize:
|
||||
m.optimize(optimizer)
|
||||
|
|
@ -316,23 +316,22 @@ def toy_ARD(max_iters=1000, kernel_type='linear', num_samples=300, D=4, optimize
|
|||
Y /= Y.std()
|
||||
|
||||
if kernel_type == 'linear':
|
||||
kernel = GPy.kern.linear(X.shape[1], ARD=1)
|
||||
kernel = GPy.kern.Linear(X.shape[1], ARD=1)
|
||||
elif kernel_type == 'rbf_inv':
|
||||
kernel = GPy.kern.rbf_inv(X.shape[1], ARD=1)
|
||||
kernel = GPy.kern.RBF_inv(X.shape[1], ARD=1)
|
||||
else:
|
||||
kernel = GPy.kern.rbf(X.shape[1], ARD=1)
|
||||
kernel += GPy.kern.white(X.shape[1]) + GPy.kern.bias(X.shape[1])
|
||||
kernel = GPy.kern.RBF(X.shape[1], ARD=1)
|
||||
kernel += GPy.kern.White(X.shape[1]) + GPy.kern.Bias(X.shape[1])
|
||||
m = GPy.models.GPRegression(X, Y, kernel)
|
||||
# len_prior = GPy.priors.inverse_gamma(1,18) # 1, 25
|
||||
# m.set_prior('.*lengthscale',len_prior)
|
||||
|
||||
if optimize:
|
||||
m.optimize(optimizer='scg', max_iters=max_iters, messages=1)
|
||||
m.optimize(optimizer='scg', max_iters=max_iters)
|
||||
|
||||
if plot:
|
||||
m.kern.plot_ARD()
|
||||
|
||||
print m
|
||||
return m
|
||||
|
||||
def toy_ARD_sparse(max_iters=1000, kernel_type='linear', num_samples=300, D=4, optimize=True, plot=True):
|
||||
|
|
@ -355,36 +354,39 @@ def toy_ARD_sparse(max_iters=1000, kernel_type='linear', num_samples=300, D=4, o
|
|||
Y /= Y.std()
|
||||
|
||||
if kernel_type == 'linear':
|
||||
kernel = GPy.kern.linear(X.shape[1], ARD=1)
|
||||
kernel = GPy.kern.Linear(X.shape[1], ARD=1)
|
||||
elif kernel_type == 'rbf_inv':
|
||||
kernel = GPy.kern.rbf_inv(X.shape[1], ARD=1)
|
||||
kernel = GPy.kern.RBF_inv(X.shape[1], ARD=1)
|
||||
else:
|
||||
kernel = GPy.kern.rbf(X.shape[1], ARD=1)
|
||||
kernel += GPy.kern.bias(X.shape[1])
|
||||
kernel = GPy.kern.RBF(X.shape[1], ARD=1)
|
||||
#kernel += GPy.kern.Bias(X.shape[1])
|
||||
X_variance = np.ones(X.shape) * 0.5
|
||||
m = GPy.models.SparseGPRegression(X, Y, kernel, X_variance=X_variance)
|
||||
# len_prior = GPy.priors.inverse_gamma(1,18) # 1, 25
|
||||
# m.set_prior('.*lengthscale',len_prior)
|
||||
|
||||
if optimize:
|
||||
m.optimize(optimizer='scg', max_iters=max_iters, messages=1)
|
||||
m.optimize(optimizer='scg', max_iters=max_iters)
|
||||
|
||||
if plot:
|
||||
m.kern.plot_ARD()
|
||||
|
||||
print m
|
||||
return m
|
||||
|
||||
def robot_wireless(max_iters=100, kernel=None, optimize=True, plot=True):
|
||||
"""Predict the location of a robot given wirelss signal strength readings."""
|
||||
data = GPy.util.datasets.robot_wireless()
|
||||
try:import pods
|
||||
except ImportError:
|
||||
print 'pods unavailable, see https://github.com/sods/ods for example datasets'
|
||||
return
|
||||
data = pods.datasets.robot_wireless()
|
||||
|
||||
# create simple GP Model
|
||||
m = GPy.models.GPRegression(data['Y'], data['X'], kernel=kernel)
|
||||
|
||||
# optimize
|
||||
if optimize:
|
||||
m.optimize(messages=True, max_iters=max_iters)
|
||||
m.optimize(max_iters=max_iters)
|
||||
|
||||
Xpredict = m.predict(data['Ytest'])[0]
|
||||
if plot:
|
||||
|
|
@ -396,13 +398,16 @@ def robot_wireless(max_iters=100, kernel=None, optimize=True, plot=True):
|
|||
|
||||
sse = ((data['Xtest'] - Xpredict)**2).sum()
|
||||
|
||||
print m
|
||||
print('Sum of squares error on test data: ' + str(sse))
|
||||
return m
|
||||
|
||||
def silhouette(max_iters=100, optimize=True, plot=True):
|
||||
"""Predict the pose of a figure given a silhouette. This is a task from Agarwal and Triggs 2004 ICML paper."""
|
||||
data = GPy.util.datasets.silhouette()
|
||||
try:import pods
|
||||
except ImportError:
|
||||
print 'pods unavailable, see https://github.com/sods/ods for example datasets'
|
||||
return
|
||||
data = pods.datasets.silhouette()
|
||||
|
||||
# create simple GP Model
|
||||
m = GPy.models.GPRegression(data['X'], data['Y'])
|
||||
|
|
@ -414,32 +419,38 @@ def silhouette(max_iters=100, optimize=True, plot=True):
|
|||
print m
|
||||
return m
|
||||
|
||||
def sparse_GP_regression_1D(num_samples=400, num_inducing=5, max_iters=100, optimize=True, plot=True):
|
||||
def sparse_GP_regression_1D(num_samples=400, num_inducing=5, max_iters=100, optimize=True, plot=True, checkgrad=False):
|
||||
"""Run a 1D example of a sparse GP regression."""
|
||||
# sample inputs and outputs
|
||||
X = np.random.uniform(-3., 3., (num_samples, 1))
|
||||
Y = np.sin(X) + np.random.randn(num_samples, 1) * 0.05
|
||||
# construct kernel
|
||||
rbf = GPy.kern.rbf(1)
|
||||
rbf = GPy.kern.RBF(1)
|
||||
# create simple GP Model
|
||||
m = GPy.models.SparseGPRegression(X, Y, kernel=rbf, num_inducing=num_inducing)
|
||||
m.checkgrad(verbose=1)
|
||||
|
||||
if checkgrad:
|
||||
m.checkgrad()
|
||||
|
||||
if optimize:
|
||||
m.optimize('tnc', messages=1, max_iters=max_iters)
|
||||
m.optimize('tnc', max_iters=max_iters)
|
||||
|
||||
if plot:
|
||||
m.plot()
|
||||
|
||||
return m
|
||||
|
||||
def sparse_GP_regression_2D(num_samples=400, num_inducing=50, max_iters=100, optimize=True, plot=True):
|
||||
def sparse_GP_regression_2D(num_samples=400, num_inducing=50, max_iters=100, optimize=True, plot=True, nan=False):
|
||||
"""Run a 2D example of a sparse GP regression."""
|
||||
np.random.seed(1234)
|
||||
X = np.random.uniform(-3., 3., (num_samples, 2))
|
||||
Y = np.sin(X[:, 0:1]) * np.sin(X[:, 1:2]) + np.random.randn(num_samples, 1) * 0.05
|
||||
if nan:
|
||||
inan = np.random.binomial(1,.2,size=Y.shape)
|
||||
Y[inan] = np.nan
|
||||
|
||||
# construct kernel
|
||||
rbf = GPy.kern.rbf(2)
|
||||
rbf = GPy.kern.RBF(2)
|
||||
|
||||
# create simple GP Model
|
||||
m = GPy.models.SparseGPRegression(X, Y, kernel=rbf, num_inducing=num_inducing)
|
||||
|
|
@ -462,7 +473,7 @@ def sparse_GP_regression_2D(num_samples=400, num_inducing=50, max_iters=100, opt
|
|||
|
||||
def uncertain_inputs_sparse_regression(max_iters=200, 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))
|
||||
fig, axes = pb.subplots(1, 2, figsize=(12, 5), sharex=True, sharey=True)
|
||||
|
||||
# sample inputs and outputs
|
||||
S = np.ones((20, 1))
|
||||
|
|
@ -471,8 +482,7 @@ def uncertain_inputs_sparse_regression(max_iters=200, optimize=True, plot=True):
|
|||
# likelihood = GPy.likelihoods.Gaussian(Y)
|
||||
Z = np.random.uniform(-3., 3., (7, 1))
|
||||
|
||||
k = GPy.kern.rbf(1)
|
||||
|
||||
k = GPy.kern.RBF(1)
|
||||
# create simple GP Model - no input uncertainty on this one
|
||||
m = GPy.models.SparseGPRegression(X, Y, kernel=k, Z=Z)
|
||||
|
||||
|
|
@ -485,7 +495,7 @@ def uncertain_inputs_sparse_regression(max_iters=200, optimize=True, plot=True):
|
|||
print m
|
||||
|
||||
# the same Model with uncertainty
|
||||
m = GPy.models.SparseGPRegression(X, Y, kernel=k, Z=Z, X_variance=S)
|
||||
m = GPy.models.SparseGPRegression(X, Y, kernel=GPy.kern.RBF(1), Z=Z, X_variance=S)
|
||||
if optimize:
|
||||
m.optimize('scg', messages=1, max_iters=max_iters)
|
||||
if plot:
|
||||
|
|
|
|||
|
|
@ -1,37 +0,0 @@
|
|||
# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
|
||||
# Licensed under the BSD 3-clause license (see LICENSE.txt)
|
||||
|
||||
import pylab as pb
|
||||
import numpy as np
|
||||
import GPy
|
||||
|
||||
def toy_1d(optimize=True, plot=True):
|
||||
N = 2000
|
||||
M = 20
|
||||
|
||||
#create data
|
||||
X = np.linspace(0,32,N)[:,None]
|
||||
Z = np.linspace(0,32,M)[:,None]
|
||||
Y = np.sin(X) + np.cos(0.3*X) + np.random.randn(*X.shape)/np.sqrt(50.)
|
||||
|
||||
m = GPy.models.SVIGPRegression(X,Y, batchsize=10, Z=Z)
|
||||
m.constrain_bounded('noise_variance',1e-3,1e-1)
|
||||
m.constrain_bounded('white_variance',1e-3,1e-1)
|
||||
|
||||
m.param_steplength = 1e-4
|
||||
|
||||
if plot:
|
||||
fig = pb.figure()
|
||||
ax = fig.add_subplot(111)
|
||||
def cb(foo):
|
||||
ax.cla()
|
||||
m.plot(ax=ax,Z_height=-3)
|
||||
ax.set_ylim(-3,3)
|
||||
fig.canvas.draw()
|
||||
|
||||
if optimize:
|
||||
m.optimize(500, callback=cb, callback_interval=1)
|
||||
|
||||
if plot:
|
||||
m.plot_traces()
|
||||
return m
|
||||
|
|
@ -1,153 +0,0 @@
|
|||
# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
|
||||
# Licensed under the BSD 3-clause license (see LICENSE.txt)
|
||||
|
||||
|
||||
"""
|
||||
Code of Tutorials
|
||||
"""
|
||||
|
||||
import pylab as pb
|
||||
pb.ion()
|
||||
import numpy as np
|
||||
import GPy
|
||||
|
||||
def tuto_GP_regression(optimize=True, plot=True):
|
||||
"""The detailed explanations of the commands used in this file can be found in the tutorial section"""
|
||||
|
||||
X = np.random.uniform(-3.,3.,(20,1))
|
||||
Y = np.sin(X) + np.random.randn(20,1)*0.05
|
||||
|
||||
kernel = GPy.kern.rbf(input_dim=1, variance=1., lengthscale=1.)
|
||||
|
||||
m = GPy.models.GPRegression(X, Y, kernel)
|
||||
|
||||
print m
|
||||
if plot:
|
||||
m.plot()
|
||||
|
||||
m.constrain_positive('')
|
||||
|
||||
m.unconstrain('') # may be used to remove the previous constrains
|
||||
m.constrain_positive('.*rbf_variance')
|
||||
m.constrain_bounded('.*lengthscale',1.,10. )
|
||||
m.constrain_fixed('.*noise',0.0025)
|
||||
|
||||
if optimize:
|
||||
m.optimize()
|
||||
m.optimize_restarts(num_restarts = 10)
|
||||
|
||||
#######################################################
|
||||
#######################################################
|
||||
# sample inputs and outputs
|
||||
X = np.random.uniform(-3.,3.,(50,2))
|
||||
Y = np.sin(X[:,0:1]) * np.sin(X[:,1:2])+np.random.randn(50,1)*0.05
|
||||
|
||||
# define kernel
|
||||
ker = GPy.kern.Matern52(2,ARD=True) + GPy.kern.white(2)
|
||||
|
||||
# create simple GP model
|
||||
m = GPy.models.GPRegression(X, Y, ker)
|
||||
|
||||
# contrain all parameters to be positive
|
||||
m.constrain_positive('')
|
||||
|
||||
# optimize and plot
|
||||
if optimize:
|
||||
m.optimize('tnc', max_f_eval = 1000)
|
||||
if plot:
|
||||
m.plot()
|
||||
|
||||
print m
|
||||
return(m)
|
||||
|
||||
def tuto_kernel_overview(optimize=True, plot=True):
|
||||
"""The detailed explanations of the commands used in this file can be found in the tutorial section"""
|
||||
ker1 = GPy.kern.rbf(1) # Equivalent to ker1 = GPy.kern.rbf(input_dim=1, variance=1., lengthscale=1.)
|
||||
ker2 = GPy.kern.rbf(input_dim=1, variance = .75, lengthscale=2.)
|
||||
ker3 = GPy.kern.rbf(1, .5, .5)
|
||||
|
||||
print ker2
|
||||
|
||||
if plot:
|
||||
ker1.plot()
|
||||
ker2.plot()
|
||||
ker3.plot()
|
||||
|
||||
k1 = GPy.kern.rbf(1,1.,2.)
|
||||
k2 = GPy.kern.Matern32(1, 0.5, 0.2)
|
||||
|
||||
# Product of kernels
|
||||
k_prod = k1.prod(k2) # By default, tensor=False
|
||||
k_prodtens = k1.prod(k2,tensor=True)
|
||||
|
||||
# Sum of kernels
|
||||
k_add = k1.add(k2) # By default, tensor=False
|
||||
k_addtens = k1.add(k2,tensor=True)
|
||||
|
||||
k1 = GPy.kern.rbf(1,1.,2)
|
||||
k2 = GPy.kern.periodic_Matern52(1,variance=1e3, lengthscale=1, period = 1.5, lower=-5., upper = 5)
|
||||
|
||||
k = k1 * k2 # equivalent to k = k1.prod(k2)
|
||||
print k
|
||||
|
||||
# Simulate sample paths
|
||||
X = np.linspace(-5,5,501)[:,None]
|
||||
Y = np.random.multivariate_normal(np.zeros(501),k.K(X),1)
|
||||
|
||||
k1 = GPy.kern.rbf(1)
|
||||
k2 = GPy.kern.Matern32(1)
|
||||
k3 = GPy.kern.white(1)
|
||||
|
||||
k = k1 + k2 + k3
|
||||
print k
|
||||
|
||||
k.constrain_positive('.*var')
|
||||
k.constrain_fixed(np.array([1]),1.75)
|
||||
k.tie_params('.*len')
|
||||
k.unconstrain('white')
|
||||
k.constrain_bounded('white',lower=1e-5,upper=.5)
|
||||
print k
|
||||
|
||||
k_cst = GPy.kern.bias(1,variance=1.)
|
||||
k_mat = GPy.kern.Matern52(1,variance=1., lengthscale=3)
|
||||
Kanova = (k_cst + k_mat).prod(k_cst + k_mat,tensor=True)
|
||||
print Kanova
|
||||
|
||||
# sample inputs and outputs
|
||||
X = np.random.uniform(-3.,3.,(40,2))
|
||||
Y = 0.5*X[:,:1] + 0.5*X[:,1:] + 2*np.sin(X[:,:1]) * np.sin(X[:,1:])
|
||||
|
||||
# Create GP regression model
|
||||
m = GPy.models.GPRegression(X, Y, Kanova)
|
||||
|
||||
if plot:
|
||||
fig = pb.figure(figsize=(5,5))
|
||||
ax = fig.add_subplot(111)
|
||||
m.plot(ax=ax)
|
||||
|
||||
pb.figure(figsize=(20,3))
|
||||
pb.subplots_adjust(wspace=0.5)
|
||||
axs = pb.subplot(1,5,1)
|
||||
m.plot(ax=axs)
|
||||
pb.subplot(1,5,2)
|
||||
pb.ylabel("= ",rotation='horizontal',fontsize='30')
|
||||
axs = pb.subplot(1,5,3)
|
||||
m.plot(ax=axs, which_parts=[False,True,False,False])
|
||||
pb.ylabel("cst +",rotation='horizontal',fontsize='30')
|
||||
axs = pb.subplot(1,5,4)
|
||||
m.plot(ax=axs, which_parts=[False,False,True,False])
|
||||
pb.ylabel("+ ",rotation='horizontal',fontsize='30')
|
||||
axs = pb.subplot(1,5,5)
|
||||
pb.ylabel("+ ",rotation='horizontal',fontsize='30')
|
||||
m.plot(ax=axs, which_parts=[False,False,False,True])
|
||||
|
||||
return(m)
|
||||
|
||||
|
||||
def model_interaction(optimize=True, plot=True):
|
||||
X = np.random.randn(20,1)
|
||||
Y = np.sin(X) + np.random.randn(*X.shape)*0.01 + 5.
|
||||
k = GPy.kern.rbf(1) + GPy.kern.bias(1)
|
||||
m = GPy.models.GPRegression(X, Y, kernel=k)
|
||||
return m
|
||||
|
||||
Loading…
Add table
Add a link
Reference in a new issue