apunkt/GPy
1
0
Fork 0
mirror of https://github.com/SheffieldML/GPy.git synced 2026-04-27 22:06:22 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions

fixed the multiple optima demo

Browse source
This commit is contained in:
James Hensman 2013-06-05 13:08:51 +01:00
parent aef3ba72ab
commit adb8a98cb3
2 changed files with 24 additions and 29 deletions
Download patch file Download diff file Expand all files Collapse all files

2
GPy/core/parameterised.py
View file

@ -85,7 +85,7 @@ class parameterised(object):
else: else:
return self._get_params()[matches] return self._get_params()[matches]
else: else:
raise AttributeError, "no parameter matches %s" % name raise AttributeError, "no parameter matches %s" % regexp
def __setitem__(self, name, val): def __setitem__(self, name, val):
""" """

51
GPy/examples/regression.py
View file

@ -181,7 +181,7 @@ def coregionalisation_sparse(optim_iters=100):
return m return m
def multiple_optima(gene_number=937,resolution=80, model_restarts=10, seed=10000, optim_iters=100): def multiple_optima(gene_number=937,resolution=80, model_restarts=10, seed=10000, optim_iters=300):
"""Show an example of a multimodal error surface for Gaussian process regression. Gene 939 has bimodal behaviour where the noisey mode is higher.""" """Show an example of a multimodal error surface for Gaussian process regression. Gene 939 has bimodal behaviour where the noisey mode is higher."""
# Contour over a range of length scales and signal/noise ratios. # Contour over a range of length scales and signal/noise ratios.
@ -197,7 +197,7 @@ def multiple_optima(gene_number=937,resolution=80, model_restarts=10, seed=10000
data['Y'] = data['Y'] - np.mean(data['Y']) 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)
pb.contour(length_scales, log_SNRs, np.exp(lls), 20) pb.contour(length_scales, log_SNRs, np.exp(lls), 20, cmap=pb.cm.jet)
ax = pb.gca() ax = pb.gca()
pb.xlabel('length scale') pb.xlabel('length scale')
pb.ylabel('log_10 SNR') pb.ylabel('log_10 SNR')
@ -211,18 +211,20 @@ def multiple_optima(gene_number=937,resolution=80, model_restarts=10, seed=10000
optim_point_y = np.empty(2) optim_point_y = np.empty(2)
np.random.seed(seed=seed) np.random.seed(seed=seed)
for i in range(0, model_restarts): for i in range(0, model_restarts):
kern = GPy.kern.rbf(1, variance=np.random.exponential(1.), lengthscale=np.random.exponential(50.)) + GPy.kern.white(1,variance=np.random.exponential(1.)) #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.GP_regression(data['X'],data['Y'], kernel=kern) m = GPy.models.GP_regression(data['X'],data['Y'], kernel=kern)
optim_point_x[0] = m.get('rbf_lengthscale') m['noise_variance'] = np.random.uniform(1e-3,1)
optim_point_y[0] = np.log10(m.get('rbf_variance')) - np.log10(m.get('white_variance')); optim_point_x[0] = m['rbf_lengthscale']
optim_point_y[0] = np.log10(m['rbf_variance']) - np.log10(m['noise_variance']);
# optimize # optimize
m.ensure_default_constraints() m.ensure_default_constraints()
m.optimize(xtol=1e-6, ftol=1e-6, max_f_eval=optim_iters) m.optimize('scg', xtol=1e-6, ftol=1e-6, max_f_eval=optim_iters)
optim_point_x[1] = m.get('rbf_lengthscale') optim_point_x[1] = m['rbf_lengthscale']
optim_point_y[1] = np.log10(m.get('rbf_variance')) - np.log10(m.get('white_variance')); optim_point_y[1] = np.log10(m['rbf_variance']) - np.log10(m['noise_variance']);
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') 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')
models.append(m) models.append(m)
@ -231,39 +233,32 @@ def multiple_optima(gene_number=937,resolution=80, model_restarts=10, seed=10000
ax.set_ylim(ylim) ax.set_ylim(ylim)
return (models, lls) return (models, lls)
def _contour_data(data, length_scales, log_SNRs, signal_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. """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. :data_set: A data set from the utils.datasets director.
:length_scales: a list of length scales to explore for the contour plot. :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. :log_SNRs: a list of base 10 logarithm signal to noise ratios to explore for the contour plot.
:signal_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 = [] lls = []
total_var = np.var(data['Y']) total_var = np.var(data['Y'])
kernel = kernel_call(1, variance=1., lengthscale=1.)
model = GPy.models.GP_regression(data['X'], data['Y'], kernel=kernel)
for log_SNR in log_SNRs: for log_SNR in log_SNRs:
SNR = 10**log_SNR SNR = 10.**log_SNR
noise_var = total_var/(1.+SNR)
signal_var = total_var - noise_var
model.kern['.*variance'] = signal_var
model['noise_variance'] = noise_var
length_scale_lls = [] length_scale_lls = []
for length_scale in length_scales: for length_scale in length_scales:
noise_var = 1. model['.*lengthscale'] = length_scale
signal_var = SNR
noise_var = noise_var/(noise_var + signal_var)*total_var
signal_var = signal_var/(noise_var + signal_var)*total_var
signal_kernel = signal_kernel_call(1, variance=signal_var, lengthscale=length_scale)
noise_kernel = GPy.kern.white(1, variance=noise_var)
kernel = signal_kernel + noise_kernel
K = kernel.K(data['X'])
total_var = (np.dot(np.dot(data['Y'].T,GPy.util.linalg.pdinv(K)[0]), data['Y'])/data['Y'].shape[0])[0,0]
noise_var *= total_var
signal_var *= total_var
kernel = signal_kernel_call(1, variance=signal_var, lengthscale=length_scale) + GPy.kern.white(1, variance=noise_var)
model = GPy.models.GP_regression(data['X'], data['Y'], kernel=kernel)
model.constrain_positive('')
length_scale_lls.append(model.log_likelihood()) length_scale_lls.append(model.log_likelihood())
lls.append(length_scale_lls) lls.append(length_scale_lls)
return np.array(lls) return np.array(lls)
def sparse_GP_regression_1D(N = 400, M = 5, optim_iters=100): def sparse_GP_regression_1D(N = 400, M = 5, optim_iters=100):