title: Multi-output spatial statistics with gaussian processes
body: <p>I've been investigating Gaussian processes lately. The perspective of probabilistic multi-output is promising in my field. In particular, spatial statistics. But I encountered three problems:</p>
  
  <ol>
  <li>multi-ouput</li>
  <li>overfitting and</li>
  <li>anisotropy.</li>
  </ol>
  
  <p>Let me run a simple case study with the <code>meuse</code> data set (from the R package <code>sp</code>).</p>
  
  <p><strong>UPDATE</strong>: The Jupyter notebook used for this question, and updated according to <a href="https://stackoverflow.com/a/43622022/1298376">Grr's answer</a>, is <a href="https://gist.github.com/essicolo/8aa8768ba24323d921f285db1e9d8272" rel="noreferrer">here</a>.</p>
  
  <pre><code>import pandas as pd
  import numpy as np
  import matplotlib.pylab as plt
  %matplotlib inline
  
  meuse = pd.read_csv(filepath_or_buffer='https://gist.githubusercontent.com/essicolo/91a2666f7c5972a91bca763daecdc5ff/raw/056bda04114d55b793469b2ab0097ec01a6d66c6/meuse.csv', sep=',')
  </code></pre>
  
  <p>For the example, we will focus on copper and lead.</p>
  
  <pre><code>fig = plt.figure(figsize=(12,8))
  ax1 = fig.add_subplot(121, aspect=1)
  ax1.set_title('Lead')
  ax1.scatter(x=meuse.x, y=meuse.y, s=meuse.lead, alpha=0.5, color='grey')
  
  ax2 = fig.add_subplot(122, aspect=1)
  ax2.set_title('Copper')
  ax2.scatter(x=meuse.x, y=meuse.y, s=meuse.copper, alpha=0.5, color='orange')
  </code></pre>
  
  <p><a href="https://i.stack.imgur.com/9aVqI.png" rel="noreferrer"><img src="https://i.stack.imgur.com/9aVqI.png" alt="enter image description here"></a></p>
  
  <p>In fact, concentrations of copper and lead are correlated.</p>
  
  <pre><code>plt.plot(meuse['lead'], meuse['copper'], '.')
  plt.xlabel('Lead')
  plt.ylabel('Copper')
  </code></pre>
  
  <p><a href="https://i.stack.imgur.com/sAjxO.png" rel="noreferrer"><img src="https://i.stack.imgur.com/sAjxO.png" alt="enter image description here"></a></p>
  
  <p>This is thus a multi-output problem.</p>
  
  <pre><code>from sklearn.gaussian_process.kernels import RBF
  from sklearn.gaussian_process import GaussianProcessRegressor as GPR
  reg = GPR(kernel=RBF())
  reg.fit(X=meuse[['x', 'y']], y=meuse[['lead', 'copper']])
  predicted = reg.predict(meuse[['x', 'y']])
  </code></pre>
  
  <p>First question: <strong>Is the kernel built for correlated multi-output when y has more than one dimension?</strong> If not, how may I specify the kernel?</p>
  
  <p>I keep on with the analysis to show the second problem, <strong>overfitting</strong>: </p>
  
  <pre><code>fig = plt.figure(figsize=(12,4))
  ax1 = fig.add_subplot(121)
  ax1.set_title('Lead')
  ax1.set_xlabel('Measured')
  ax1.set_ylabel('Predicted')
  ax1.plot(meuse.lead, predicted[:,0], '.')
  
  ax2 = fig.add_subplot(122)
  ax2.set_title('Copper')
  ax2.set_xlabel('Measured')
  ax2.set_ylabel('Predicted')
  ax2.plot(meuse.copper, predicted[:,1], '.')
  </code></pre>
  
  <p><a href="https://i.stack.imgur.com/W2GjM.png" rel="noreferrer"><img src="https://i.stack.imgur.com/W2GjM.png" alt="enter image description here"></a></p>
  
  <p>I created a grid of x and y coordinates and all concentrations on that grid were predicted as zeros.</p>
  
  <p>Finally, a last concern which particularly arises in 3D for soils: how may I specify <strong>anisotropy</strong> in such models?</p>