nc, likfunc, x, y, z);
  f = [m+2*sqrt(s2); flipdim(m-2*sqrt(s2),1)];
  fill([z; flipdim(z,1)], f, [7 7 7]/8)
  hold on; plot(z, m); plot(x, y, '+')
</pre>

<p>First, we guess that a <a href="../cov/covSEiso.m">squared
exponential</a> covariance function <tt>covSEiso</tt> may be suitable.
This covariance function takes two hyperparameters: a characteristic
length-scale and a signal standard deviation (magnitude). These
hyperparameters are non-negative and represented by their logarithms;
thus, initializing 