regression.html

高斯过程在回归和分类问题中的应用

  1.5461
  0.1443
</pre>

Note that the hyperparameters learned here a close, but not identical to the
log hyperparameters <tt>1.0, 1.0, 0.1</tt> used when generating the data. The
discrepancy is partially due to the small training sample size, and partially
due to the fact that we only get information about the process in a very
limited range of input values. Predicting with the learned hyperparameters
produces:<br>

<center><img src="figlf.gif"></center><br>

which is not very different from the plot above (based on the hyperparameters
used to generate the data). Repeating the experiment with more training points
distributed over a wider range leads to more accurate estimates.</p>

Note, that above we have use the same functional form for the covariance
function, as was used to generate the data. In practise things are seldom so
simple, and one may have to try different covariance functions. Here, we try
to explore how a Matern form, with a shape parameter of 3/2 
(see eq. 4.17) would do.

<pre>
  covfunc2 = {'covSum',{'covMatern3iso','covNoise'}};
  loghyper2 = minimize([-1; -1; -1], 'gpr', -100, covfunc2, x, y);
</pre>

Comparing the value of the marginal likelihoods for the two models

<pre>
  -gpr(loghyper, covfunc, x, y)
  -gpr(loghyper2, covfunc2, x, y)
</pre>

with values of <tt>-15.6</tt> for SE and <tt>-18.0</tt> for Matern3, shows that
the SE covariance function is about <tt>exp(18.0-15.6)=11</tt> times more
probable than the Matern form for these data (in agreement with the data
generating process). The predictions from the worse Matern-based model

<pre>
  [mu S2] = gpr(loghyper2, covfunc2, x, y, xstar);
  S2 = S2 - exp(2*loghyper2(3));
</pre>  

look like this:<br>

<center><img src="figlm.gif"></center><br>

Notice how the uncertainty grows more rapidly in the vicinity of data-points,
reflecting the property that for the Matern class 
with a shape parameter of 3/2 the stochastic process is not twice mean-square
differentiable (and is thus
much less smooth than the SE covariance function).


<h3 id="ard">3. GPR demo using gpr.m with multi-dimensional input space and
ARD</h3>

This demonstration illustrates the use of a Gaussian Process
regression for a multi-dimensional input, and illustrates
the use of automatic relevance determination (ARD). It is based
on <a href="#williams-rasmussen-96">Williams and Rasmussen (1996)</a>.</p>

You can either follow the example below or run the script <a
href="../gpml-demo/demo_gparm.m">demo_gparm.m</a>.</p>

We initially consider a 2-d nonlinear robot arm mapping problem as defined in
<a href="#mackay-92">Mackay (1992)</a>

<pre>
  f(x_1,x_2) = r_1 cos (x_1) + r_2 cos(x_1 + x_2)
</pre>

where x_1 is chosen randomly in [-1.932, -0.453], x_2 is chosen randomly in
[0.534, 3.142], and r_1 = 2.0, r_2 = 1.3.  The target values are obtained by
adding Gaussian noise of variance 0.0025 to f(x_1,x_2).  Following <a
href="#neal-96">Neal (1996)</a> we add four further inputs, two of which (x_3
and x_4) are copies of x_1 and x_2 corrupted by additive Gaussian noise of
standard deviation 0.02, and two of which (x_5 and x_6) are N(0,1) Gaussian
noise variables.  Our dataset has n=200 training points and nstar=200 test
points.</p>

The training and test data is contained in the file <a
href="../gpml-demo/data_6darm.mat">data_6darm.mat</a>. The raw training data is
in the input matrix <tt>X</tt> (<tt>n=</tt>200 by <tt>D=</tt>6) and the target
vector <tt>y</tt> (200 by 1). Assuming that the current directory is
<tt>gpml-demo</tt> we need to add the path of the code, and load the data:

<pre>
  addpath('../gpml');    % add path for the code
  load data_6darm
</pre>

We first check the scaling of these variables using

<pre>
  mean(X), std(X), mean(y), std(y)
</pre>

In particular we might be concerned if the standard deviation is very different
for different input dimensions; however, that is not the case here so we do not
carry out rescaling for <tt>X</tt>.  However, <tt>y</tt> has a non-zero mean
which is not appropriate if we assume a zero-mean GP. We could add a constant
onto the SE covariance function corresponding to a prior on constant offsets,
but here instead we centre <tt>y</tt> by setting:

<pre>
  offset = mean(y);
  y = y - offset;         % centre targets around 0.
</pre>

We use Gaussian process regression with a squared exponential covariance
function, and allow a separate lengthscale for each input dimension, as in
eqs. 5.1 and 5.2.  These lengthscales (and the other hyperparameters sigma_f
and sigma_n) are adapted by maximizing the marginal likelihood (eq. 5.8)
w.r.t. the hyperparameters. The covariance function is specified by

<pre>
  covfunc = {'covSum', {'covSEard','covNoise'}};
</pre>

We now wish to train the GP by optimizing the hyperparameters.  The
hyperparameters are stored as <tt>logtheta = [log(ell_1), log(ell_2),
... log(ell_6), log(sigma_f), log(sigma_n)]</tt> (as <tt>D = 6</tt>), and are
initialized to <tt>logtheta0 = [0; 0; 0; 0; 0; 0; 0; log(sqrt(0.1))]</tt>.  The
last value means that the initial noise variance sigma^2_n is set to 0.1.  The
marginal likelihood is optimized using the <tt>minimize</tt> function.

<pre>
  logtheta0 = [0; 0; 0; 0; 0; 0; 0; log(sqrt(0.1))];
  [logtheta, fvals, iter] = minimize(logtheta0, 'gpr', -100, covfunc, X, y);
  exp(logtheta) 
</pre>

By doing <tt>exp(logtheta)</tt> we find that:

<pre>
  ell_1      1.804377
  ell_2      1.963956
  ell_3      8.884361
  ell_4     34.417657
  ell_5   1081.610451
  ell_6    375.445823
  sigma_f    2.379139
  sigma_n    0.050835
</pre>

Notice that the lengthscales ell_1 and ell_2 are short indicating that inputs
x_1 and x_2 are relevant to the task.  The noisy inputs x_3 and x_4 have longer
lengthscales, indicating they are less relevant, and the pure noise inputs x_5
and x_6 have very long lengthscales, so they are effectively irrelevant to the
problem, as indeed we would hope.  The process std deviation sigma_f is similar
in magnitude to the standard deviation of the data <tt>std(y) = 1.2186</tt>.
The learned noise standard deviation sigma_n is almost exactly equal to the
generative noise level sqrt(0.0025)=0.05.</p>

We now make predictions on the test points and assess the accuracy of these
predictions. The test data is in <tt>Xstar</tt> and <tt>ystar</tt>. Recall that
the training targets were centered, thus we must adjust the predictions by
<tt>offset</tt>:

<pre>
  [fstar S2] = gpr(logtheta, covfunc, X, y, Xstar);
  fstar = fstar + offset;  % add back offset to get true prediction
</pre>

We compute the residuals, the mean squared error (mse) and the predictive log
likelihood (pll). Note that the predictive variance for <tt>ystar</tt> includes
the noise variance, as explained on p. 18.

<pre>
  res = fstar-ystar;  % residuals
  mse = mean(res.^2)
  pll = -0.5*mean(log(2*pi*S2)+res.^2./S2)
</pre>

The mean squared error is 0.002489. Note that this is almost equal to the value
0.0025, as would be obtained from the perfect predictor, due to the added noise
with variance 0.0025.</p>
 
We can also plot the residuals and the predictive (noise-free) variance for
each test case. Note that the order along the x-axis is arbitrary.

<pre>
  subplot(211), plot(res,'.'), ylabel('residuals'), xlabel('test case')
  subplot(212), plot(sqrt(S2),'.'), ylabel('predictive std deviation'), xlabel('test case')
</pre>


<h3 id="ref">4. References</h3>

<ul>
<li id=mackay-92>D. J. C. MacKay, <a
href="http://www.inference.phy.cam.ac.uk/mackay/backprop.nc.ps.gz">A Practical
Bayesian Framework for Backpropagation Networks</a>, Neural Computation 4,
448-472, (1992)
</li>
<li id=neal-96>
R. M. Neal, <a href="http://www.cs.utoronto.ca/~radford/bnn.book.html">Bayesian
Learning for Neural Networks</a>, Springer, (1996)
</li>
<li id=williams-rasmussen-96>
C. K. I. Williams and C. E. Rasmussen, <a
href="http://www.gaussianprocess.org/#williams-rasmussen-96">Gaussian Processes
for Regression</a>, Advances in Neural Information Processing Systems 8, pp
514-520, MIT Press (1996).
</li>
</ul>

Go back to the <a href="http://www.gaussianprocess.org/gpml">web page</a> for
Gaussian Processes for Machine Learning.

<hr>
<!-- Created: Mon Oct 17 14:16:13 CEST 2005 -->
<!-- hhmts start -->
Last modified: Mon Mar 27 16:10:14 CEST 2006
<!-- hhmts end -->
  </body>
</html>