[WIP] Gaussian Process revisited (refactoring of correlation models + some new features/bugfixes)

This PR contains a proposal for a refactoring of the gaussian_process subpackage (mainly the correlation_models module) and some novel features which become possible because of this refactoring. Moreover, some subtle bugs are fixed. In general, I tried to keep the external interface of GaussianProcess fixed and to change only internals.The main purpose of this PR is to begin a discussion on the future of the GP subpackage and revive its further development.

CHANGES:
 * Correlation models are now classes instead of functions and inherit from abstract base class StationaryCorrelation. This reduces code-redundancy (DRY) and separates GPs and their correlation models more strictly. Moreover, non-stationary correlation models become possible (see examples/gaussian_process/plot_gp_nonstationary.py) This was not possible formerly since only the differences of the datapoints were passed to the correlation models but not the datapoints themselves.
 * The hyperparameters theta are estimated as maximum-a-posterior estimates rather than maximum-likelihood estimates. At the moment, the prior on the hyperparameters theta is uniform such that both result in the same value. However, a non-uniform prior could be used in the future in cases where the risk of overfitting theta to the data is serious (e.g., forcing the factor analysis distance to be close to diagonal etc.)
 * Hyperparameters theta are represented as 1D ndarrays instead of 2D 1xn arrays in GaussianProcess. 

NEW FEATURES:
 * Matern correlation models for nu=1.5 and nu=2.5 have been added (see https://en.wikipedia.org/wiki/Mat%C3%A9rn_covariance_function). An example script showing the potential benefit of the Matern correlation model compared to squared-exponential and absolute-exponential was added under examples/gaussian_process/plot_matern_kernel.py (see attached image)
 * squared_exponential, absolute_exponential and Matern correlation models support factor analysis distance. This can be seen as an extension of learning dimension-specific length scales in which also correlations of feature dimensions can be taken into account. See Rasmussen and Williams 2006, p107 for details. An example script showing the potential benefit of this extension was added under examples/gaussian_process/plot_gp_learning_curve.py (see attached image). 
 * GP supports additional optimizers for ML estimation passed as callables. In addition to 'fmin_cobyla' and 'Welch', GaussianProcess allows now to pass other optimizers directly as callables via the parameter optimizer.

BUGFIXES
* Setting optimal_rlf_value correctly in GPs and actually maximizing likelihood over several random starts. The sign of the optimal_rlf_value was set wrongly in _arg_max_reduced_likelihood_function() since it was set to the negative of the return of reduced_likelihood_function(). This error was probably caused by confusing reduced_likelihood_function(theta)  with minus_reduced_likelihood_function(log10t). It resulted, however, in _arg_max_reduced_likelihood_function() returning the worst local maxima instead of the best one when performing several random_starts.
* A typo in gp_diabetes_dataset.py was fixed which caused a crash of the example (gp.theta_ versus gp.theta)

To be discussed:
 * The factor analysis distance has hyperparameters that can can take on arbitrary real values (not just positive ones). Since sklearn enforces the hyperparameters theta to be positive, this is internally handled by taking the log of the corresponding components of theta. Are there better ideas?


![plot_gp_learning_curve](https://f.cloud.github.com/assets/1116263/2310367/a7499dde-a2e2-11e3-80f5-6e98a23f0ee8.png)
![plot_matern_kernel](https://f.cloud.github.com/assets/1116263/2310372/adf004fc-a2e2-11e3-8d25-85e69e09fd11.png)


Is there any way we could avoid overriding the __call__ method, and call transform() explicitly instead by setting properties on calls to fit()?