Show Sidebar Hide Sidebar

# Gaussian Rrocess Regression with Noise-Level Estimation in Scikit-learn

This example illustrates that GPR with a sum-kernel including a WhiteKernel can estimate the noise level of data. An illustration of the log-marginal-likelihood (LML) landscape shows that there exist two local maxima of LML. The first corresponds to a model with a high noise level and a large length scale, which explains all variations in the data by noise. The second one has a smaller noise level and shorter length scale, which explains most of the variation by the noise-free functional relationship. The second model has a higher likelihood; however, depending on the initial value for the hyperparameters, the gradient-based optimization might also converge to the high-noise solution. It is thus important to repeat the optimization several times for different initializations.

#### New to Plotly?¶

You can set up Plotly to work in online or offline mode, or in jupyter notebooks.
We also have a quick-reference cheatsheet (new!) to help you get started!

### Version¶

In [1]:
import sklearn
sklearn.__version__

Out[1]:
'0.18.1'

### Imports¶

This tutorial imports GaussianProcessRegressor, RBF and WhiteKernel.

In [2]:
import plotly.plotly as py
import plotly.graph_objs as go

import numpy as np
from matplotlib import pyplot as plt
from matplotlib.colors import LogNorm

from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, WhiteKernel


### Calculations¶

In [3]:
rng = np.random.RandomState(0)
X = rng.uniform(0, 5, 20)[:, np.newaxis]
y = 0.5 * np.sin(3 * X[:, 0]) + rng.normal(0, 0.5, X.shape[0])


### First Run¶

In [4]:
kernel = 1.0 * RBF(length_scale=100.0, length_scale_bounds=(1e-2, 1e3)) \
+ WhiteKernel(noise_level=1, noise_level_bounds=(1e-10, 1e+1))
gp = GaussianProcessRegressor(kernel=kernel,
alpha=0.0).fit(X, y)
X_ = np.linspace(0, 5, 100)
y_mean, y_cov = gp.predict(X_[:, np.newaxis], return_cov=True)

trace1 = go.Scatter(x=X_,y=y_mean,
mode='lines',
line=dict(color='black'),
showlegend=False,
fill = 'tonexty'
)

trace2 = go.Scatter(x=X_, y=y_mean - np.sqrt(np.diag(y_cov)),
mode='lines',
line=dict(color='black'),
showlegend=False,
fill='tonexty'
)

trace3 = go.Scatter(x=X_, y=y_mean + np.sqrt(np.diag(y_cov)),
mode='lines',
line=dict(color='black'),
showlegend=False
)

trace4 = go.Scatter(x=X_, y=0.5*np.sin(3*X_),
mode='lines',
line=dict(color='red'),
showlegend=False
)

trace5 = go.Scatter(x=X[:, 0], y=y,
mode='markers',
marker=dict(color='red'),
showlegend=False
)
data = [ trace2, trace3,  trace1,trace4, trace5]
layout = go.Layout(title="Initial: %s<br>Optimum: %s<br>Log-Marginal-Likelihood: %s"
% (kernel, gp.kernel_,
gp.log_marginal_likelihood(gp.kernel_.theta)))
fig = go.Figure(data=data, layout=layout)

In [5]:
py.iplot(fig)

Out[5]:

### Second Run¶

In [6]:
kernel = 1.0 * RBF(length_scale=1.0, length_scale_bounds=(1e-2, 1e3)) \
+ WhiteKernel(noise_level=1e-5, noise_level_bounds=(1e-10, 1e+1))
gp = GaussianProcessRegressor(kernel=kernel,
alpha=0.0).fit(X, y)
X_ = np.linspace(0, 5, 100)
y_mean, y_cov = gp.predict(X_[:, np.newaxis], return_cov=True)

trace1 = go.Scatter(x=X_,y=y_mean,
mode='lines',
line=dict(color='black'),
showlegend=False,
fill = 'tozeroy'
)

trace2 = go.Scatter(x=X_, y=y_mean - np.sqrt(np.diag(y_cov)),
mode='lines',
line=dict(color='black'),
showlegend=False,
fill='tonexty'
)

trace3 = go.Scatter(x=X_, y=y_mean + np.sqrt(np.diag(y_cov)),
mode='lines',
line=dict(color='black'),
showlegend=False,
fill='tonexty'
)

trace4 = go.Scatter(x=X_, y=0.5*np.sin(3*X_),
mode='lines',
line=dict(color='red'),
showlegend=False
)

trace5 = go.Scatter(x=X[:, 0], y=y,
mode='markers',
marker=dict(color='red'),
showlegend=False
)
data = [ trace2, trace3, trace1, trace4, trace5]
layout = go.Layout(title="Initial: %s<br>Optimum: %s<br>Log-Marginal-Likelihood: %s"
% (kernel, gp.kernel_,
gp.log_marginal_likelihood(gp.kernel_.theta)))
fig = go.Figure(data=data, layout=layout)

In [7]:
py.iplot(fig)

Out[7]:

### Plot LML landscape¶

In [8]:
theta0 = np.logspace(-2, 3, 49)
theta1 = np.logspace(-2, 0, 50)
Theta0, Theta1 = np.meshgrid(theta0, theta1)
LML = [[gp.log_marginal_likelihood(np.log([0.36, Theta0[i, j], Theta1[i, j]]))
for i in range(Theta0.shape[0])] for j in range(Theta0.shape[1])]
LML = np.array(LML).T

vmin, vmax = (-LML).min(), (-LML).max()
vmax = 50

trace = go.Contour(x=theta0,
y=theta1,
z=-LML,
ncontours=np.logspace(np.log10(vmin), np.log10(vmax), 50),
contours=dict(coloring='lines')
)
layout = go.Layout(title="Log-marginal-likelihood",
xaxis=dict(type='log', title="Length-scale",
showgrid=False),
yaxis=dict(type='log', title="Noise-level",
showgrid=False)
)
fig = go.Figure(data=[trace], layout=layout)

In [10]:
py.iplot(fig)

Out[10]:

Authors:

    Jan Hendrik Metzen <jhm@informatik.uni-bremen.de>



    BSD 3 clause
Still need help?