diff --git a/GPy/core/parameterised.py b/GPy/core/parameterised.py index c8d8ce4d..e091e820 100644 --- a/GPy/core/parameterised.py +++ b/GPy/core/parameterised.py @@ -85,7 +85,7 @@ class parameterised(object): else: return self._get_params()[matches] else: - raise AttributeError, "no parameter matches %s" % name + raise AttributeError, "no parameter matches %s" % regexp def __setitem__(self, name, val): """ diff --git a/GPy/examples/regression.py b/GPy/examples/regression.py index 310a0691..7989b3c1 100644 --- a/GPy/examples/regression.py +++ b/GPy/examples/regression.py @@ -181,7 +181,7 @@ def coregionalisation_sparse(optim_iters=100): 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.""" # 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']) 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() pb.xlabel('length scale') 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) 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.)) + 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) - optim_point_x[0] = m.get('rbf_lengthscale') - optim_point_y[0] = np.log10(m.get('rbf_variance')) - np.log10(m.get('white_variance')); + 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']); # optimize 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_y[1] = np.log10(m.get('rbf_variance')) - np.log10(m.get('white_variance')); + optim_point_x[1] = m['rbf_lengthscale'] + 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') models.append(m) @@ -231,39 +233,32 @@ def multiple_optima(gene_number=937,resolution=80, model_restarts=10, seed=10000 ax.set_ylim(ylim) 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. :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. - :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 = [] 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: - 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 = [] + for length_scale in length_scales: - noise_var = 1. - 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('') + model['.*lengthscale'] = length_scale length_scale_lls.append(model.log_likelihood()) + lls.append(length_scale_lls) + return np.array(lls) def sparse_GP_regression_1D(N = 400, M = 5, optim_iters=100):