def multiple_optima(gene_number=<NUM_LIT>, resolution=<NUM_LIT>, model_restarts=<NUM_LIT:10>, seed=<NUM_LIT>, max_iters=<NUM_LIT>, optimize=True, plot=True):
<EOL>length_scales = np.linspace(<NUM_LIT:0.1>, <NUM_LIT>, resolution)<EOL>log_SNRs = np.linspace(-<NUM_LIT>, <NUM_LIT>, resolution)<EOL>try:import pods<EOL>except ImportError:<EOL><INDENT>print('<STR_LIT>')<EOL>return<EOL><DEDENT>data = pods.datasets.della_gatta_TRP63_gene_expression(data_set='<STR_LIT>',gene_number=gene_number)<EOL>data['<STR_LIT:Y>'] = data['<STR_LIT:Y>'] - np.mean(data['<STR_LIT:Y>'])<EOL>lls = GPy.examples.regression._contour_data(data, length_scales, log_SNRs, GPy.kern.RBF)<EOL>if plot:<EOL><INDENT>pb.contour(length_scales, log_SNRs, np.exp(lls), <NUM_LIT:20>, cmap=pb.cm.jet)<EOL>ax = pb.gca()<EOL>pb.xlabel('<STR_LIT>')<EOL>pb.ylabel('<STR_LIT>')<EOL>xlim = ax.get_xlim()<EOL>ylim = ax.get_ylim()<EOL><DEDENT>models = []<EOL>optim_point_x = np.empty(<NUM_LIT:2>)<EOL>optim_point_y = np.empty(<NUM_LIT:2>)<EOL>np.random.seed(seed=seed)<EOL>for i in range(<NUM_LIT:0>, model_restarts):<EOL><INDENT>kern = GPy.kern.RBF(<NUM_LIT:1>, variance=np.random.uniform(<NUM_LIT>, <NUM_LIT:1>), lengthscale=np.random.uniform(<NUM_LIT:5>, <NUM_LIT:50>))<EOL>m = GPy.models.GPRegression(data['<STR_LIT:X>'], data['<STR_LIT:Y>'], kernel=kern)<EOL>m.likelihood.variance = np.random.uniform(<NUM_LIT>, <NUM_LIT:1>)<EOL>optim_point_x[<NUM_LIT:0>] = m.rbf.lengthscale<EOL>optim_point_y[<NUM_LIT:0>] = np.log10(m.rbf.variance) - np.log10(m.likelihood.variance);<EOL>if optimize:<EOL><INDENT>m.optimize('<STR_LIT>', xtol=<NUM_LIT>, ftol=<NUM_LIT>, max_iters=max_iters)<EOL><DEDENT>optim_point_x[<NUM_LIT:1>] = m.rbf.lengthscale<EOL>optim_point_y[<NUM_LIT:1>] = np.log10(m.rbf.variance) - np.log10(m.likelihood.variance);<EOL>if plot:<EOL><INDENT>pb.arrow(optim_point_x[<NUM_LIT:0>], optim_point_y[<NUM_LIT:0>], optim_point_x[<NUM_LIT:1>] - optim_point_x[<NUM_LIT:0>], optim_point_y[<NUM_LIT:1>] - optim_point_y[<NUM_LIT:0>], label=str(i), head_length=<NUM_LIT:1>, head_width=<NUM_LIT:0.5>, fc='<STR_LIT:k>', ec='<STR_LIT:k>')<EOL><DEDENT>models.append(m)<EOL><DEDENT>if plot:<EOL><INDENT>ax.set_xlim(xlim)<EOL>ax.set_ylim(ylim)<EOL><DEDENT>return m<EOL>
Show an example of a multimodal error surface for Gaussian process
regression. Gene 939 has bimodal behaviour where the noisy mode is
higher.