1# Copyright (c) 2012-2014, GPy authors (see AUTHORS.txt). 2# Licensed under the BSD 3-clause license (see LICENSE.txt) 3 4""" 5Gaussian Processes regression examples 6""" 7try: 8 from matplotlib import pyplot as pb 9except: 10 pass 11import numpy as np 12import GPy 13 14def olympic_marathon_men(optimize=True, plot=True): 15 """Run a standard Gaussian process regression on the Olympic marathon data.""" 16 try:import pods 17 except ImportError: 18 print('pods unavailable, see https://github.com/sods/ods for example datasets') 19 return 20 data = pods.datasets.olympic_marathon_men() 21 22 # create simple GP Model 23 m = GPy.models.GPRegression(data['X'], data['Y']) 24 25 # set the lengthscale to be something sensible (defaults to 1) 26 m.kern.lengthscale = 10. 27 28 if optimize: 29 m.optimize('bfgs', max_iters=200) 30 if plot: 31 m.plot(plot_limits=(1850, 2050)) 32 33 return m 34 35def coregionalization_toy(optimize=True, plot=True): 36 """ 37 A simple demonstration of coregionalization on two sinusoidal functions. 38 """ 39 #build a design matrix with a column of integers indicating the output 40 X1 = np.random.rand(50, 1) * 8 41 X2 = np.random.rand(30, 1) * 5 42 43 #build a suitable set of observed variables 44 Y1 = np.sin(X1) + np.random.randn(*X1.shape) * 0.05 45 Y2 = np.sin(X2) + np.random.randn(*X2.shape) * 0.05 + 2. 46 47 m = GPy.models.GPCoregionalizedRegression(X_list=[X1,X2], Y_list=[Y1,Y2]) 48 49 if optimize: 50 m.optimize('bfgs', max_iters=100) 51 52 if plot: 53 slices = GPy.util.multioutput.get_slices([X1,X2]) 54 m.plot(fixed_inputs=[(1,0)],which_data_rows=slices[0],Y_metadata={'output_index':0}) 55 m.plot(fixed_inputs=[(1,1)],which_data_rows=slices[1],Y_metadata={'output_index':1},ax=pb.gca()) 56 return m 57 58def coregionalization_sparse(optimize=True, plot=True): 59 """ 60 A simple demonstration of coregionalization on two sinusoidal functions using sparse approximations. 61 """ 62 #build a design matrix with a column of integers indicating the output 63 X1 = np.random.rand(50, 1) * 8 64 X2 = np.random.rand(30, 1) * 5 65 66 #build a suitable set of observed variables 67 Y1 = np.sin(X1) + np.random.randn(*X1.shape) * 0.05 68 Y2 = np.sin(X2) + np.random.randn(*X2.shape) * 0.05 + 2. 69 70 m = GPy.models.SparseGPCoregionalizedRegression(X_list=[X1,X2], Y_list=[Y1,Y2]) 71 72 if optimize: 73 m.optimize('bfgs', max_iters=100) 74 75 if plot: 76 slices = GPy.util.multioutput.get_slices([X1,X2]) 77 m.plot(fixed_inputs=[(1,0)],which_data_rows=slices[0],Y_metadata={'output_index':0}) 78 m.plot(fixed_inputs=[(1,1)],which_data_rows=slices[1],Y_metadata={'output_index':1},ax=pb.gca()) 79 pb.ylim(-3,) 80 81 return m 82 83def epomeo_gpx(max_iters=200, optimize=True, plot=True): 84 """ 85 Perform Gaussian process regression on the latitude and longitude data 86 from the Mount Epomeo runs. Requires gpxpy to be installed on your system 87 to load in the data. 88 """ 89 try:import pods 90 except ImportError: 91 print('pods unavailable, see https://github.com/sods/ods for example datasets') 92 return 93 data = pods.datasets.epomeo_gpx() 94 num_data_list = [] 95 for Xpart in data['X']: 96 num_data_list.append(Xpart.shape[0]) 97 98 num_data_array = np.array(num_data_list) 99 num_data = num_data_array.sum() 100 Y = np.zeros((num_data, 2)) 101 t = np.zeros((num_data, 2)) 102 start = 0 103 for Xpart, index in zip(data['X'], range(len(data['X']))): 104 end = start+Xpart.shape[0] 105 t[start:end, :] = np.hstack((Xpart[:, 0:1], 106 index*np.ones((Xpart.shape[0], 1)))) 107 Y[start:end, :] = Xpart[:, 1:3] 108 109 num_inducing = 200 110 Z = np.hstack((np.linspace(t[:,0].min(), t[:, 0].max(), num_inducing)[:, None], 111 np.random.randint(0, 4, num_inducing)[:, None])) 112 113 k1 = GPy.kern.RBF(1) 114 k2 = GPy.kern.Coregionalize(output_dim=5, rank=5) 115 k = k1**k2 116 117 m = GPy.models.SparseGPRegression(t, Y, kernel=k, Z=Z, normalize_Y=True) 118 m.constrain_fixed('.*variance', 1.) 119 m.inducing_inputs.constrain_fixed() 120 m.Gaussian_noise.variance.constrain_bounded(1e-3, 1e-1) 121 m.optimize(max_iters=max_iters,messages=True) 122 123 return m 124 125def multiple_optima(gene_number=937, resolution=80, model_restarts=10, seed=10000, max_iters=300, optimize=True, plot=True): 126 """ 127 Show an example of a multimodal error surface for Gaussian process 128 regression. Gene 939 has bimodal behaviour where the noisy mode is 129 higher. 130 """ 131 132 # Contour over a range of length scales and signal/noise ratios. 133 length_scales = np.linspace(0.1, 60., resolution) 134 log_SNRs = np.linspace(-3., 4., resolution) 135 136 try:import pods 137 except ImportError: 138 print('pods unavailable, see https://github.com/sods/ods for example datasets') 139 return 140 data = pods.datasets.della_gatta_TRP63_gene_expression(data_set='della_gatta',gene_number=gene_number) 141 # data['Y'] = data['Y'][0::2, :] 142 # data['X'] = data['X'][0::2, :] 143 144 data['Y'] = data['Y'] - np.mean(data['Y']) 145 146 lls = GPy.examples.regression._contour_data(data, length_scales, log_SNRs, GPy.kern.RBF) 147 if plot: 148 pb.contour(length_scales, log_SNRs, np.exp(lls), 20, cmap=pb.cm.jet) 149 ax = pb.gca() 150 pb.xlabel('length scale') 151 pb.ylabel('log_10 SNR') 152 153 xlim = ax.get_xlim() 154 ylim = ax.get_ylim() 155 156 # Now run a few optimizations 157 models = [] 158 optim_point_x = np.empty(2) 159 optim_point_y = np.empty(2) 160 np.random.seed(seed=seed) 161 for i in range(0, model_restarts): 162 # kern = GPy.kern.RBF(1, variance=np.random.exponential(1.), lengthscale=np.random.exponential(50.)) 163 kern = GPy.kern.RBF(1, variance=np.random.uniform(1e-3, 1), lengthscale=np.random.uniform(5, 50)) 164 165 m = GPy.models.GPRegression(data['X'], data['Y'], kernel=kern) 166 m.likelihood.variance = np.random.uniform(1e-3, 1) 167 optim_point_x[0] = m.rbf.lengthscale 168 optim_point_y[0] = np.log10(m.rbf.variance) - np.log10(m.likelihood.variance); 169 170 # optimize 171 if optimize: 172 m.optimize('scg', xtol=1e-6, ftol=1e-6, max_iters=max_iters) 173 174 optim_point_x[1] = m.rbf.lengthscale 175 optim_point_y[1] = np.log10(m.rbf.variance) - np.log10(m.likelihood.variance); 176 177 if plot: 178 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') 179 models.append(m) 180 181 if plot: 182 ax.set_xlim(xlim) 183 ax.set_ylim(ylim) 184 return m # (models, lls) 185 186def _contour_data(data, length_scales, log_SNRs, kernel_call=GPy.kern.RBF): 187 """ 188 Evaluate the GP objective function for a given data set for a range of 189 signal to noise ratios and a range of lengthscales. 190 191 :data_set: A data set from the utils.datasets director. 192 :length_scales: a list of length scales to explore for the contour plot. 193 :log_SNRs: a list of base 10 logarithm signal to noise ratios to explore for the contour plot. 194 :kernel: a kernel to use for the 'signal' portion of the data. 195 """ 196 197 lls = [] 198 total_var = np.var(data['Y']) 199 kernel = kernel_call(1, variance=1., lengthscale=1.) 200 model = GPy.models.GPRegression(data['X'], data['Y'], kernel=kernel) 201 for log_SNR in log_SNRs: 202 SNR = 10.**log_SNR 203 noise_var = total_var / (1. + SNR) 204 signal_var = total_var - noise_var 205 model.kern['.*variance'] = signal_var 206 model.likelihood.variance = noise_var 207 length_scale_lls = [] 208 209 for length_scale in length_scales: 210 model['.*lengthscale'] = length_scale 211 length_scale_lls.append(model.log_likelihood()) 212 213 lls.append(length_scale_lls) 214 215 return np.array(lls) 216 217 218def olympic_100m_men(optimize=True, plot=True): 219 """Run a standard Gaussian process regression on the Rogers and Girolami olympics data.""" 220 try:import pods 221 except ImportError: 222 print('pods unavailable, see https://github.com/sods/ods for example datasets') 223 return 224 data = pods.datasets.olympic_100m_men() 225 226 # create simple GP Model 227 m = GPy.models.GPRegression(data['X'], data['Y']) 228 229 # set the lengthscale to be something sensible (defaults to 1) 230 m.rbf.lengthscale = 10 231 232 if optimize: 233 m.optimize('bfgs', max_iters=200) 234 235 if plot: 236 m.plot(plot_limits=(1850, 2050)) 237 return m 238 239def toy_rbf_1d(optimize=True, plot=True): 240 """Run a simple demonstration of a standard Gaussian process fitting it to data sampled from an RBF covariance.""" 241 try:import pods 242 except ImportError: 243 print('pods unavailable, see https://github.com/sods/ods for example datasets') 244 return 245 data = pods.datasets.toy_rbf_1d() 246 247 # create simple GP Model 248 m = GPy.models.GPRegression(data['X'], data['Y']) 249 250 if optimize: 251 m.optimize('bfgs') 252 if plot: 253 m.plot() 254 255 return m 256 257def toy_rbf_1d_50(optimize=True, plot=True): 258 """Run a simple demonstration of a standard Gaussian process fitting it to data sampled from an RBF covariance.""" 259 try:import pods 260 except ImportError: 261 print('pods unavailable, see https://github.com/sods/ods for example datasets') 262 return 263 data = pods.datasets.toy_rbf_1d_50() 264 265 # create simple GP Model 266 m = GPy.models.GPRegression(data['X'], data['Y']) 267 268 if optimize: 269 m.optimize('bfgs') 270 if plot: 271 m.plot() 272 273 return m 274 275def toy_poisson_rbf_1d_laplace(optimize=True, plot=True): 276 """Run a simple demonstration of a standard Gaussian process fitting it to data sampled from an RBF covariance.""" 277 optimizer='scg' 278 x_len = 100 279 X = np.linspace(0, 10, x_len)[:, None] 280 f_true = np.random.multivariate_normal(np.zeros(x_len), GPy.kern.RBF(1).K(X)) 281 Y = np.array([np.random.poisson(np.exp(f)) for f in f_true])[:,None] 282 283 kern = GPy.kern.RBF(1) 284 poisson_lik = GPy.likelihoods.Poisson() 285 laplace_inf = GPy.inference.latent_function_inference.Laplace() 286 287 # create simple GP Model 288 m = GPy.core.GP(X, Y, kernel=kern, likelihood=poisson_lik, inference_method=laplace_inf) 289 290 if optimize: 291 m.optimize(optimizer) 292 if plot: 293 m.plot() 294 # plot the real underlying rate function 295 pb.plot(X, np.exp(f_true), '--k', linewidth=2) 296 297 return m 298 299def toy_ARD(max_iters=1000, kernel_type='linear', num_samples=300, D=4, optimize=True, plot=True): 300 # Create an artificial dataset where the values in the targets (Y) 301 # only depend in dimensions 1 and 3 of the inputs (X). Run ARD to 302 # see if this dependency can be recovered 303 X1 = np.sin(np.sort(np.random.rand(num_samples, 1) * 10, 0)) 304 X2 = np.cos(np.sort(np.random.rand(num_samples, 1) * 10, 0)) 305 X3 = np.exp(np.sort(np.random.rand(num_samples, 1), 0)) 306 X4 = np.log(np.sort(np.random.rand(num_samples, 1), 0)) 307 X = np.hstack((X1, X2, X3, X4)) 308 309 Y1 = np.asarray(2 * X[:, 0] + 3).reshape(-1, 1) 310 Y2 = np.asarray(4 * (X[:, 2] - 1.5 * X[:, 0])).reshape(-1, 1) 311 Y = np.hstack((Y1, Y2)) 312 313 Y = np.dot(Y, np.random.rand(2, D)); 314 Y = Y + 0.2 * np.random.randn(Y.shape[0], Y.shape[1]) 315 Y -= Y.mean() 316 Y /= Y.std() 317 318 if kernel_type == 'linear': 319 kernel = GPy.kern.Linear(X.shape[1], ARD=1) 320 elif kernel_type == 'rbf_inv': 321 kernel = GPy.kern.RBF_inv(X.shape[1], ARD=1) 322 else: 323 kernel = GPy.kern.RBF(X.shape[1], ARD=1) 324 kernel += GPy.kern.White(X.shape[1]) + GPy.kern.Bias(X.shape[1]) 325 m = GPy.models.GPRegression(X, Y, kernel) 326 # len_prior = GPy.priors.inverse_gamma(1,18) # 1, 25 327 # m.set_prior('.*lengthscale',len_prior) 328 329 if optimize: 330 m.optimize(optimizer='scg', max_iters=max_iters) 331 332 if plot: 333 m.kern.plot_ARD() 334 335 return m 336 337def toy_ARD_sparse(max_iters=1000, kernel_type='linear', num_samples=300, D=4, optimize=True, plot=True): 338 # Create an artificial dataset where the values in the targets (Y) 339 # only depend in dimensions 1 and 3 of the inputs (X). Run ARD to 340 # see if this dependency can be recovered 341 X1 = np.sin(np.sort(np.random.rand(num_samples, 1) * 10, 0)) 342 X2 = np.cos(np.sort(np.random.rand(num_samples, 1) * 10, 0)) 343 X3 = np.exp(np.sort(np.random.rand(num_samples, 1), 0)) 344 X4 = np.log(np.sort(np.random.rand(num_samples, 1), 0)) 345 X = np.hstack((X1, X2, X3, X4)) 346 347 Y1 = np.asarray(2 * X[:, 0] + 3)[:, None] 348 Y2 = np.asarray(4 * (X[:, 2] - 1.5 * X[:, 0]))[:, None] 349 Y = np.hstack((Y1, Y2)) 350 351 Y = np.dot(Y, np.random.rand(2, D)); 352 Y = Y + 0.2 * np.random.randn(Y.shape[0], Y.shape[1]) 353 Y -= Y.mean() 354 Y /= Y.std() 355 356 if kernel_type == 'linear': 357 kernel = GPy.kern.Linear(X.shape[1], ARD=1) 358 elif kernel_type == 'rbf_inv': 359 kernel = GPy.kern.RBF_inv(X.shape[1], ARD=1) 360 else: 361 kernel = GPy.kern.RBF(X.shape[1], ARD=1) 362 #kernel += GPy.kern.Bias(X.shape[1]) 363 X_variance = np.ones(X.shape) * 0.5 364 m = GPy.models.SparseGPRegression(X, Y, kernel, X_variance=X_variance) 365 # len_prior = GPy.priors.inverse_gamma(1,18) # 1, 25 366 # m.set_prior('.*lengthscale',len_prior) 367 368 if optimize: 369 m.optimize(optimizer='scg', max_iters=max_iters) 370 371 if plot: 372 m.kern.plot_ARD() 373 374 return m 375 376def robot_wireless(max_iters=100, kernel=None, optimize=True, plot=True): 377 """Predict the location of a robot given wirelss signal strength readings.""" 378 try:import pods 379 except ImportError: 380 print('pods unavailable, see https://github.com/sods/ods for example datasets') 381 return 382 data = pods.datasets.robot_wireless() 383 384 # create simple GP Model 385 m = GPy.models.GPRegression(data['Y'], data['X'], kernel=kernel) 386 387 # optimize 388 if optimize: 389 m.optimize(max_iters=max_iters) 390 391 Xpredict = m.predict(data['Ytest'])[0] 392 if plot: 393 pb.plot(data['Xtest'][:, 0], data['Xtest'][:, 1], 'r-') 394 pb.plot(Xpredict[:, 0], Xpredict[:, 1], 'b-') 395 pb.axis('equal') 396 pb.title('WiFi Localization with Gaussian Processes') 397 pb.legend(('True Location', 'Predicted Location')) 398 399 sse = ((data['Xtest'] - Xpredict)**2).sum() 400 401 print(('Sum of squares error on test data: ' + str(sse))) 402 return m 403 404def silhouette(max_iters=100, optimize=True, plot=True): 405 """Predict the pose of a figure given a silhouette. This is a task from Agarwal and Triggs 2004 ICML paper.""" 406 try:import pods 407 except ImportError: 408 print('pods unavailable, see https://github.com/sods/ods for example datasets') 409 return 410 data = pods.datasets.silhouette() 411 412 # create simple GP Model 413 m = GPy.models.GPRegression(data['X'], data['Y']) 414 415 # optimize 416 if optimize: 417 m.optimize(messages=True, max_iters=max_iters) 418 419 print(m) 420 return m 421 422def sparse_GP_regression_1D(num_samples=400, num_inducing=5, max_iters=100, optimize=True, plot=True, checkgrad=False): 423 """Run a 1D example of a sparse GP regression.""" 424 # sample inputs and outputs 425 X = np.random.uniform(-3., 3., (num_samples, 1)) 426 Y = np.sin(X) + np.random.randn(num_samples, 1) * 0.05 427 # construct kernel 428 rbf = GPy.kern.RBF(1) 429 # create simple GP Model 430 m = GPy.models.SparseGPRegression(X, Y, kernel=rbf, num_inducing=num_inducing) 431 432 if checkgrad: 433 m.checkgrad() 434 435 if optimize: 436 m.optimize('tnc', max_iters=max_iters) 437 438 if plot: 439 m.plot() 440 441 return m 442 443def sparse_GP_regression_2D(num_samples=400, num_inducing=50, max_iters=100, optimize=True, plot=True, nan=False): 444 """Run a 2D example of a sparse GP regression.""" 445 np.random.seed(1234) 446 X = np.random.uniform(-3., 3., (num_samples, 2)) 447 Y = np.sin(X[:, 0:1]) * np.sin(X[:, 1:2]) + np.random.randn(num_samples, 1) * 0.05 448 if nan: 449 inan = np.random.binomial(1,.2,size=Y.shape) 450 Y[inan] = np.nan 451 452 # construct kernel 453 rbf = GPy.kern.RBF(2) 454 455 # create simple GP Model 456 m = GPy.models.SparseGPRegression(X, Y, kernel=rbf, num_inducing=num_inducing) 457 458 # contrain all parameters to be positive (but not inducing inputs) 459 m['.*len'] = 2. 460 461 m.checkgrad() 462 463 # optimize 464 if optimize: 465 m.optimize('tnc', messages=1, max_iters=max_iters) 466 467 # plot 468 if plot: 469 m.plot() 470 471 print(m) 472 return m 473 474def uncertain_inputs_sparse_regression(max_iters=200, optimize=True, plot=True): 475 """Run a 1D example of a sparse GP regression with uncertain inputs.""" 476 fig, axes = pb.subplots(1, 2, figsize=(12, 5), sharex=True, sharey=True) 477 478 # sample inputs and outputs 479 S = np.ones((20, 1)) 480 X = np.random.uniform(-3., 3., (20, 1)) 481 Y = np.sin(X) + np.random.randn(20, 1) * 0.05 482 # likelihood = GPy.likelihoods.Gaussian(Y) 483 Z = np.random.uniform(-3., 3., (7, 1)) 484 485 k = GPy.kern.RBF(1) 486 # create simple GP Model - no input uncertainty on this one 487 m = GPy.models.SparseGPRegression(X, Y, kernel=k, Z=Z) 488 489 if optimize: 490 m.optimize('scg', messages=1, max_iters=max_iters) 491 492 if plot: 493 m.plot(ax=axes[0]) 494 axes[0].set_title('no input uncertainty') 495 print(m) 496 497 # the same Model with uncertainty 498 m = GPy.models.SparseGPRegression(X, Y, kernel=GPy.kern.RBF(1), Z=Z, X_variance=S) 499 if optimize: 500 m.optimize('scg', messages=1, max_iters=max_iters) 501 if plot: 502 m.plot(ax=axes[1]) 503 axes[1].set_title('with input uncertainty') 504 fig.canvas.draw() 505 506 print(m) 507 return m 508 509def simple_mean_function(max_iters=100, optimize=True, plot=True): 510 """ 511 The simplest possible mean function. No parameters, just a simple Sinusoid. 512 """ 513 #create simple mean function 514 mf = GPy.core.Mapping(1,1) 515 mf.f = np.sin 516 mf.update_gradients = lambda a,b: None 517 518 X = np.linspace(0,10,50).reshape(-1,1) 519 Y = np.sin(X) + 0.5*np.cos(3*X) + 0.1*np.random.randn(*X.shape) 520 521 k =GPy.kern.RBF(1) 522 lik = GPy.likelihoods.Gaussian() 523 m = GPy.core.GP(X, Y, kernel=k, likelihood=lik, mean_function=mf) 524 if optimize: 525 m.optimize(max_iters=max_iters) 526 if plot: 527 m.plot(plot_limits=(-10,15)) 528 return m 529 530def parametric_mean_function(max_iters=100, optimize=True, plot=True): 531 """ 532 A linear mean function with parameters that we'll learn alongside the kernel 533 """ 534 #create simple mean function 535 mf = GPy.core.Mapping(1,1) 536 mf.f = np.sin 537 538 X = np.linspace(0,10,50).reshape(-1,1) 539 Y = np.sin(X) + 0.5*np.cos(3*X) + 0.1*np.random.randn(*X.shape) + 3*X 540 541 mf = GPy.mappings.Linear(1,1) 542 543 k =GPy.kern.RBF(1) 544 lik = GPy.likelihoods.Gaussian() 545 m = GPy.core.GP(X, Y, kernel=k, likelihood=lik, mean_function=mf) 546 if optimize: 547 m.optimize(max_iters=max_iters) 548 if plot: 549 m.plot() 550 return m 551 552 553def warped_gp_cubic_sine(max_iters=100): 554 """ 555 A test replicating the cubic sine regression problem from 556 Snelson's paper. 557 """ 558 X = (2 * np.pi) * np.random.random(151) - np.pi 559 Y = np.sin(X) + np.random.normal(0,0.2,151) 560 Y = np.array([np.power(abs(y),float(1)/3) * (1,-1)[y<0] for y in Y]) 561 X = X[:, None] 562 Y = Y[:, None] 563 564 warp_k = GPy.kern.RBF(1) 565 warp_f = GPy.util.warping_functions.TanhFunction(n_terms=2) 566 warp_m = GPy.models.WarpedGP(X, Y, kernel=warp_k, warping_function=warp_f) 567 warp_m['.*\.d'].constrain_fixed(1.0) 568 m = GPy.models.GPRegression(X, Y) 569 m.optimize_restarts(parallel=False, robust=True, num_restarts=5, max_iters=max_iters) 570 warp_m.optimize_restarts(parallel=False, robust=True, num_restarts=5, max_iters=max_iters) 571 #m.optimize(max_iters=max_iters) 572 #warp_m.optimize(max_iters=max_iters) 573 574 print(warp_m) 575 print(warp_m['.*warp.*']) 576 577 warp_m.predict_in_warped_space = False 578 warp_m.plot(title="Warped GP - Latent space") 579 warp_m.predict_in_warped_space = True 580 warp_m.plot(title="Warped GP - Warped space") 581 m.plot(title="Standard GP") 582 warp_m.plot_warping() 583 pb.show() 584 return warp_m 585 586 587 588def multioutput_gp_with_derivative_observations(): 589 def plot_gp_vs_real(m, x, yreal, size_inputs, title, fixed_input=1, xlim=[0,11], ylim=[-1.5,3]): 590 fig, ax = pb.subplots() 591 ax.set_title(title) 592 pb.plot(x, yreal, "r", label='Real function') 593 rows = slice(0, size_inputs[0]) if fixed_input == 0 else slice(size_inputs[0], size_inputs[0]+size_inputs[1]) 594 m.plot(fixed_inputs=[(1, fixed_input)], which_data_rows=rows, xlim=xlim, ylim=ylim, ax=ax) 595 f = lambda x: np.sin(x)+0.1*(x-2.)**2-0.005*x**3 596 fd = lambda x: np.cos(x)+0.2*(x-2.)-0.015*x**2 597 N=10 # Number of observations 598 M=10 # Number of derivative observations 599 Npred=100 # Number of prediction points 600 sigma = 0.05 # Noise of observations 601 sigma_der = 0.05 # Noise of derivative observations 602 x = np.array([np.linspace(1,10,N)]).T 603 y = f(x) + np.array(sigma*np.random.normal(0,1,(N,1))) 604 605 xd = np.array([np.linspace(2,8,M)]).T 606 yd = fd(xd) + np.array(sigma_der*np.random.normal(0,1,(M,1))) 607 608 xpred = np.array([np.linspace(0,11,Npred)]).T 609 ypred_true = f(xpred) 610 ydpred_true = fd(xpred) 611 612 # squared exponential kernel: 613 se = GPy.kern.RBF(input_dim = 1, lengthscale=1.5, variance=0.2) 614 # We need to generate separate kernel for the derivative observations and give the created kernel as an input: 615 se_der = GPy.kern.DiffKern(se, 0) 616 617 #Then 618 gauss = GPy.likelihoods.Gaussian(variance=sigma**2) 619 gauss_der = GPy.likelihoods.Gaussian(variance=sigma_der**2) 620 621 # Then create the model, we give everything in lists, the order of the inputs indicates the order of the outputs 622 # Now we have the regular observations first and derivative observations second, meaning that the kernels and 623 # the likelihoods must follow the same order. Crosscovariances are automatically taken care of 624 m = GPy.models.MultioutputGP(X_list=[x, xd], Y_list=[y, yd], 625 kernel_list=[se, se_der], 626 likelihood_list=[gauss, gauss_der]) 627 628 # Optimize the model 629 m.optimize(messages=0, ipython_notebook=False) 630 631 #Plot the model, the syntax is same as for multioutput models: 632 plot_gp_vs_real(m, xpred, ydpred_true, [x.shape[0], xd.shape[0]], title='Latent function derivatives', fixed_input=1, xlim=[0,11], ylim=[-1.5,3]) 633 plot_gp_vs_real(m, xpred, ypred_true, [x.shape[0], xd.shape[0]], title='Latent function', fixed_input=0, xlim=[0,11], ylim=[-1.5,3]) 634 635 #making predictions for the values: 636 mu, var = m.predict_noiseless(Xnew=[xpred, np.empty((0,1))]) 637 638 return m 639