title: How to correctly use scikit-learn's Gaussian Process for a 2D-inputs, 1D-output regression?
body: <p>Prior to posting I did a lot of searches and found <a href="https://stackoverflow.com/questions/21320964/how-to-make-a-2d-gaussian-process-using-gpml-matlab-for-regression">this question</a> which might be exactly my problem. However, I tried what is proposed in the answer but unfortunately this did not fix it, and I couldn't add a comment to request further explanation, as I am a new member here. </p>
  
  <p>Anyway, I want to use the Gaussian Processes with scikit-learn in Python on a simple but real case to start (using the examples provided in scikit-learn's documentation). I have a 2D input set (8 couples of 2 parameters) called <strong>X</strong>. I have 8 corresponding outputs, gathered in the 1D-array <strong>y</strong>.</p>
  
  <pre><code>#  Inputs: 8 points 
  X = np.array([[p1, q1],[p2, q2],[p3, q3],[p4, q4],[p5, q5],[p6, q6],[p7, q7],[p8, q8]])
  
  # Observations: 8 couples
  y = np.array([r1,r2,r3,r4,r5,r6,r7,r8])
  </code></pre>
  
  <p>I defined an input test space <strong>x</strong>:</p>
  
  <pre><code># Input space
  x1 = np.linspace(x1min, x1max) #p
  x2 = np.linspace(x2min, x2max) #q
  x = (np.array([x1, x2])).T
  </code></pre>
  
  <p>Then I instantiate the GP model, fit it to my training data <strong>(X,y)</strong>, and make the 1D prediction <strong>y_pred</strong> on my input space <strong>x</strong>: </p>
  
  <pre><code>from sklearn.gaussian_process import GaussianProcessRegressor
  from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C
  
  kernel = C(1.0, (1e-3, 1e3)) * RBF([5,5], (1e-2, 1e2))
  gp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=15)
  gp.fit(X, y)
  y_pred, MSE = gp.predict(x, return_std=True)
  </code></pre>
  
  <p>And then I make a 3D plot: </p>
  
  <pre><code>fig = pl.figure()
  ax = fig.add_subplot(111, projection='3d')
  Xp, Yp = np.meshgrid(x1, x2)
  Zp = np.reshape(y_pred,50)
  
  surf = ax.plot_surface(Xp, Yp, Zp, rstride=1, cstride=1, cmap=cm.jet,
  linewidth=0, antialiased=False)
  pl.show()
  </code></pre>
  
  <p>This is what I obtain: </p>
  
  <p><a href="https://i.stack.imgur.com/0I6Eb.jpg" rel="noreferrer"><img src="https://i.stack.imgur.com/0I6Eb.jpg" alt="RBF[5,5]"></a></p>
  
  <p>When I modify the kernel parameters I get something like this, similar to what the poster I mentioned above got: </p>
  
  <p><a href="https://i.stack.imgur.com/3q7vn.jpg" rel="noreferrer"><img src="https://i.stack.imgur.com/3q7vn.jpg" alt="RBF[10,10]"></a></p>
  
  <p>These plots don't even match the observation from the original training points (the lower response is obtained for [65.1,37] and the highest for [92.3,54]).</p>
  
  <p>I am fairly new to GPs in 2D (also started Python not long ago) so I think I'm missing something here... Any answer would be helpful and greatly appreciated, thanks!</p>