In [1]:
%pylab inline
In [2]:
%load http://www.astroml.org/_downloads/fig_gp_example.py
In []:
"""
Gaussian Process Example
------------------------
Figure 8.10
An example of Gaussian process regression. The upper-left panel shows three
functions drawn from an unconstrained Gaussian process with squared-exponential
covari- ance of bandwidth h = 1.0. The upper-right panel adds two constraints,
and shows the 2-sigma contours of the constrained function space. The
lower-left panel shows the function space constrained by the points with error
bars. The lower-right panel shows the function space constrained by 20 noisy
points drawn from f(x) = cos(x).
"""
# Author: Jake VanderPlas
# License: BSD
# The figure produced by this code is published in the textbook
# "Statistics, Data Mining, and Machine Learning in Astronomy" (2013)
# For more information, see http://astroML.github.com
# To report a bug or issue, use the following forum:
# https://groups.google.com/forum/#!forum/astroml-general
import numpy as np
from matplotlib import pyplot as plt
from sklearn.gaussian_process import GaussianProcess
#----------------------------------------------------------------------
# This function adjusts matplotlib settings for a uniform feel in the textbook.
# Note that with usetex=True, fonts are rendered with LaTeX. This may
# result in an error if LaTeX is not installed on your system. In that case,
# you can set usetex to False.
from astroML.plotting import setup_text_plots
setup_text_plots(fontsize=8, usetex=True)
#------------------------------------------------------------
# define a squared exponential covariance function
def squared_exponential(x1, x2, h):
return np.exp(-0.5 * (x1 - x2) ** 2 / h ** 2)
#------------------------------------------------------------
# draw samples from the unconstrained covariance
np.random.seed(1)
x = np.linspace(0, 10, 100)
h = 1.0
mu = np.zeros(len(x))
C = squared_exponential(x, x[:, None], h)
draws = np.random.multivariate_normal(mu, C, 3)
#------------------------------------------------------------
# Constrain the mean and covariance with two points
x1 = np.array([2.5, 7])
y1 = np.cos(x1)
gp1 = GaussianProcess(corr='squared_exponential', theta0=0.5,
random_state=0)
gp1.fit(x1[:, None], y1)
f1, MSE1 = gp1.predict(x[:, None], eval_MSE=True)
f1_err = np.sqrt(MSE1)
#------------------------------------------------------------
# Constrain the mean and covariance with two noisy points
# scikit-learn gaussian process uses nomenclature from the geophysics
# community, where a "nugget" can be specified. The diagonal of the
# assumed covariance matrix is multiplied by the nugget. This is
# how the error on inputs is incorporated into the calculation
dy2 = 0.2
gp2 = GaussianProcess(corr='squared_exponential', theta0=0.5,
nugget=(dy2 / y1) ** 2, random_state=0)
gp2.fit(x1[:, None], y1)
f2, MSE2 = gp2.predict(x[:, None], eval_MSE=True)
f2_err = np.sqrt(MSE2)
#------------------------------------------------------------
# Constrain the mean and covariance with many noisy points
x3 = np.linspace(0, 10, 20)
y3 = np.cos(x3)
dy3 = 0.2
y3 = np.random.normal(y3, dy3)
gp3 = GaussianProcess(corr='squared_exponential', theta0=0.5,
thetaL=0.01, thetaU=10.0,
nugget=(dy3 / y3) ** 2,
random_state=0)
gp3.fit(x3[:, None], y3)
f3, MSE3 = gp3.predict(x[:, None], eval_MSE=True)
f3_err = np.sqrt(MSE3)
# we have fit for the `h` parameter: print the result here:
print "best-fit theta =", gp3.theta_[0, 0]
#------------------------------------------------------------
# Plot the diagrams
fig = plt.figure(figsize=(5, 5))
# first: plot a selection of unconstrained functions
ax = fig.add_subplot(221)
ax.plot(x, draws.T, '-k')
ax.set_ylabel('$f(x)$')
# second: plot a constrained function
ax = fig.add_subplot(222)
ax.plot(x, f1, '-', color='gray')
ax.fill_between(x, f1 - 2 * f1_err, f1 + 2 * f1_err, color='gray', alpha=0.3)
ax.plot(x1, y1, '.k', ms=6)
# third: plot a constrained function with errors
ax = fig.add_subplot(223)
ax.plot(x, f2, '-', color='gray')
ax.fill_between(x, f2 - 2 * f2_err, f2 + 2 * f2_err, color='gray', alpha=0.3)
ax.errorbar(x1, y1, dy2, fmt='.k', ms=6)
ax.set_xlabel('$x$')
ax.set_ylabel('$f(x)$')
# third: plot a more constrained function with errors
ax = fig.add_subplot(224)
ax.plot(x, f3, '-', color='gray')
ax.fill_between(x, f3 - 2 * f3_err, f3 + 2 * f3_err, color='gray', alpha=0.3)
ax.errorbar(x3, y3, dy3, fmt='.k', ms=6)
ax.plot(x, np.cos(x), ':k')
ax.set_xlabel('$x$')
for ax in fig.axes:
ax.set_xlim(0, 10)
plt.show()
In [3]:
"""
Gaussian Process Example
------------------------
Figure 8.10
An example of Gaussian process regression. The upper-left panel shows three
functions drawn from an unconstrained Gaussian process with squared-exponential
covari- ance of bandwidth h = 1.0. The upper-right panel adds two constraints,
and shows the 2-sigma contours of the constrained function space. The
lower-left panel shows the function space constrained by the points with error
bars. The lower-right panel shows the function space constrained by 20 noisy
points drawn from f(x) = cos(x).
"""
# Author: Jake VanderPlas
# License: BSD
# The figure produced by this code is published in the textbook
# "Statistics, Data Mining, and Machine Learning in Astronomy" (2013)
# For more information, see http://astroML.github.com
# To report a bug or issue, use the following forum:
# https://groups.google.com/forum/#!forum/astroml-general
import numpy as np
from matplotlib import pyplot as plt
from sklearn.gaussian_process import GaussianProcess
#----------------------------------------------------------------------
# This function adjusts matplotlib settings for a uniform feel in the textbook.
# Note that with usetex=True, fonts are rendered with LaTeX. This may
# result in an error if LaTeX is not installed on your system. In that case,
# you can set usetex to False.
from astroML.plotting import setup_text_plots
setup_text_plots(fontsize=8, usetex=True)
#------------------------------------------------------------
# define a squared exponential covariance function
def squared_exponential(x1, x2, h):
return np.exp(-0.5 * (x1 - x2) ** 2 / h ** 2)
#------------------------------------------------------------
# draw samples from the unconstrained covariance
np.random.seed(1)
x = np.linspace(0, 10, 100)
h = 1.0
mu = np.zeros(len(x))
C = squared_exponential(x, x[:, None], h)
draws = np.random.multivariate_normal(mu, C, 3)
#------------------------------------------------------------
# Constrain the mean and covariance with two points
x1 = np.array([2.5, 7])
y1 = np.cos(x1)
gp1 = GaussianProcess(corr='squared_exponential', theta0=0.5,
random_state=0)
gp1.fit(x1[:, None], y1)
f1, MSE1 = gp1.predict(x[:, None], eval_MSE=True)
f1_err = np.sqrt(MSE1)
#------------------------------------------------------------
# Constrain the mean and covariance with two noisy points
# scikit-learn gaussian process uses nomenclature from the geophysics
# community, where a "nugget" can be specified. The diagonal of the
# assumed covariance matrix is multiplied by the nugget. This is
# how the error on inputs is incorporated into the calculation
dy2 = 0.2
gp2 = GaussianProcess(corr='squared_exponential', theta0=0.5,
nugget=(dy2 / y1) ** 2, random_state=0)
gp2.fit(x1[:, None], y1)
f2, MSE2 = gp2.predict(x[:, None], eval_MSE=True)
f2_err = np.sqrt(MSE2)
#------------------------------------------------------------
# Constrain the mean and covariance with many noisy points
x3 = np.linspace(0, 10, 20)
y3 = np.cos(x3)
dy3 = 0.2
y3 = np.random.normal(y3, dy3)
gp3 = GaussianProcess(corr='squared_exponential', theta0=0.5,
thetaL=0.01, thetaU=10.0,
nugget=(dy3 / y3) ** 2,
random_state=0)
gp3.fit(x3[:, None], y3)
f3, MSE3 = gp3.predict(x[:, None], eval_MSE=True)
f3_err = np.sqrt(MSE3)
# we have fit for the `h` parameter: print the result here:
print "best-fit theta =", gp3.theta_[0, 0]
#------------------------------------------------------------
# Plot the diagrams
fig = plt.figure(figsize=(5, 5))
# first: plot a selection of unconstrained functions
ax = fig.add_subplot(221)
ax.plot(x, draws.T, '-k')
ax.set_ylabel('$f(x)$')
# second: plot a constrained function
ax = fig.add_subplot(222)
ax.plot(x, f1, '-', color='gray')
ax.fill_between(x, f1 - 2 * f1_err, f1 + 2 * f1_err, color='gray', alpha=0.3)
ax.plot(x1, y1, '.k', ms=6)
# third: plot a constrained function with errors
ax = fig.add_subplot(223)
ax.plot(x, f2, '-', color='gray')
ax.fill_between(x, f2 - 2 * f2_err, f2 + 2 * f2_err, color='gray', alpha=0.3)
ax.errorbar(x1, y1, dy2, fmt='.k', ms=6)
ax.set_xlabel('$x$')
ax.set_ylabel('$f(x)$')
# third: plot a more constrained function with errors
ax = fig.add_subplot(224)
ax.plot(x, f3, '-', color='gray')
ax.fill_between(x, f3 - 2 * f3_err, f3 + 2 * f3_err, color='gray', alpha=0.3)
ax.errorbar(x3, y3, dy3, fmt='.k', ms=6)
ax.plot(x, np.cos(x), ':k')
ax.set_xlabel('$x$')
for ax in fig.axes:
ax.set_xlim(0, 10)
plt.show()
In [4]:
import jinja2
import json
import mpld3
from mpld3 import plugins, utils
class HighlightLines(plugins.PluginBase):
"""A plugin to highlight lines on hover"""
JAVASCRIPT = """
mpld3.register_plugin("linehighlight", LineHighlightPlugin);
LineHighlightPlugin.prototype = Object.create(mpld3.Plugin.prototype);
LineHighlightPlugin.prototype.constructor = LineHighlightPlugin;
LineHighlightPlugin.prototype.requiredProps = ["line_ids"];
LineHighlightPlugin.prototype.defaultProps = {alpha_bg:0.03, alpha_fg:1.0}
function LineHighlightPlugin(fig, props){
mpld3.Plugin.call(this, fig, props);
};
LineHighlightPlugin.prototype.draw = function(){
for(var i=0; i<this.props.line_ids.length; i++){
var obj = mpld3.get_element(this.props.line_ids[i], this.fig),
alpha_fg = this.props.alpha_fg;
alpha_bg = this.props.alpha_bg;
obj.elements()
.on("mouseover", function(d, i){
d3.select(this).transition().duration(50)
.style("stroke-opacity", alpha_fg); })
.on("mouseout", function(d, i){
d3.select(this).transition().duration(200)
.style("stroke-opacity", alpha_bg); });
}
};
"""
def __init__(self, lines):
self.lines = lines
self.dict_ = {"type": "linehighlight",
"line_ids": [utils.get_id(line) for line in lines],
"alpha_bg": lines[0].get_alpha(),
"alpha_fg": 1.0}
In [9]:
"""
Enhanced Gaussian Process Example
---------------------------------
Figure 8.10
An example of Gaussian process regression. The upper-left panel shows three
functions drawn from an unconstrained Gaussian process with squared-exponential
covari- ance of bandwidth h = 1.0. The upper-right panel adds two constraints,
and shows the 2-sigma contours of the constrained function space. The
lower-left panel shows the function space constrained by the points with error
bars. The lower-right panel shows the function space constrained by 20 noisy
points drawn from f(x) = cos(x).
"""
# Original Author: Jake VanderPlas
# MPL3D mashed in by: Beth Reid & Phil Marshall at #AstroHackWeek (2014)
# License: BSD
# The figure produced by this code is published in the textbook
# "Statistics, Data Mining, and Machine Learning in Astronomy" (2013)
# For more information, see http://astroML.github.com
# To report a bug or issue, use the following forum:
# https://groups.google.com/forum/#!forum/astroml-general
import numpy as np
from matplotlib import pyplot as plt
from sklearn.gaussian_process import GaussianProcess
#----------------------------------------------------------------------
# This function adjusts matplotlib settings for a uniform feel in the textbook.
# Note that with usetex=True, fonts are rendered with LaTeX. This may
# result in an error if LaTeX is not installed on your system. In that case,
# you can set usetex to False.
from astroML.plotting import setup_text_plots
setup_text_plots(fontsize=8, usetex=True)
#------------------------------------------------------------
# define a squared exponential covariance function
def squared_exponential(x1, x2, theta0):
return np.exp(-theta0*(x1 - x2) ** 2)
#------------------------------------------------------------
# draw samples from the unconstrained covariance
np.random.seed(1)
x = np.linspace(0, 10, 100)
h = 1.0
theta0=0.5/h**2
mu = np.zeros(len(x))
C = squared_exponential(x, x[:,None], theta0)
Ndraws = 67
draws = np.random.multivariate_normal(mu, C, Ndraws)
#------------------------------------------------------------
# Constrain the mean and covariance with two points
x1 = np.array([2.5, 7])
y1 = np.cos(x1)
#y1 = 0.001*np.cos(x1)
# to get the Gaussian process to align with our definition of squared_exponential, we had to express theta0 in terms of h.
gp1 = GaussianProcess(corr='squared_exponential', theta0=theta0,
random_state=0)
gp1.fit(x1[:, None], y1)
f1, MSE1 = gp1.predict(x[:, None], eval_MSE=True)
f1_err = np.sqrt(MSE1)
K = np.zeros([len(x)+len(x1),len(x)+len(x1)])
K11 = squared_exponential(x1, x1[:,None], theta0)
K12 = squared_exponential(x, x1[:,None], theta0)
K22 = C.copy()
K[:2,:2] = K11
K[2:,2:] = K22
K[2:,:2] = K12.T
K[:2,2:] = K12
#original guess wrong?
mu1 = np.dot(K12.T,np.dot(np.linalg.inv(K11),y1))
C1 = K22 - np.dot(K12.T,np.dot(np.linalg.inv(K11),K12))
f1draws = np.random.multivariate_normal(mu1, C1, Ndraws)
In [10]:
print gp1.theta_, 0.5/h**2
print theta0
print gp1.y_std, y1.std()
In [11]:
#------------------------------------------------------------
# Constrain the mean and covariance with two noisy points
# scikit-learn gaussian process uses nomenclature from the geophysics
# community, where a "nugget" can be specified. The diagonal of the
# assumed covariance matrix is multiplied by the nugget. This is
# how the error on inputs is incorporated into the calculation
dy2 = 0.2
# documentation says added, not multiplied, so this code is wrong (?)
#original plot.
#gp2 = GaussianProcess(corr='squared_exponential', theta0=theta0,
# nugget=(dy2 / y1) ** 2, random_state=0)
#is this correct?
gp2 = GaussianProcess(corr='squared_exponential', theta0=theta0,nugget=(dy2) ** 2, random_state=0)
gp2.fit(x1[:, None], y1)
f2, MSE2 = gp2.predict(x[:, None], eval_MSE=True)
f2_err = np.sqrt(MSE2)
## now we want to make draws with this
K = np.zeros([len(x)+len(x1),len(x)+len(x1)])
K11 = squared_exponential(x1, x1[:,None], theta0) + np.eye(len(x1)) * dy2**2
K12 = squared_exponential(x, x1[:,None], theta0)
K22 = C.copy()
K[:2,:2] = K11
K[2:,2:] = K22
K[2:,:2] = K12.T
K[:2,2:] = K12
#original guess wrong?
mu2 = np.dot(K12.T,np.dot(np.linalg.inv(K11),y1))
C2 = K22 - np.dot(K12.T,np.dot(np.linalg.inv(K11),K12))
f2draws = np.random.multivariate_normal(mu2, C2, Ndraws)
#------------------------------------------------------------
# Constrain the mean and covariance with many noisy points
x3 = np.linspace(0, 10, 20)
y3 = np.cos(x3)
dy3 = 0.2
y3 = np.random.normal(y3, dy3)
gp3 = GaussianProcess(corr='squared_exponential', theta0=theta0,
thetaL=0.01, thetaU=10.0,
nugget=(dy3 / y3) ** 2,
random_state=0)
gp3.fit(x3[:, None], y3)
f3, MSE3 = gp3.predict(x[:, None], eval_MSE=True)
f3_err = np.sqrt(MSE3)
# we have fit for the `h` parameter: print the result here:
print "best-fit theta =", gp3.theta_[0, 0]
Now, we are going to animate the plots, so we need a class to enable this. This came from the http://mpld3.github.io/_downloads/random_walk.py example.
In [12]:
#------------------------------------------------------------
# Plot the diagrams
fig = plt.figure(figsize=(5, 5))
# first: plot a selection of unconstrained functions
ax = fig.add_subplot(221)
lines = ax.plot(x, draws.T, lw=4, alpha=0.03+0.3/Ndraws, color='blue')
plugins.connect(fig, HighlightLines(lines))
ax.set_ylabel('$f(x)$')
# second: plot a constrained function
ax = fig.add_subplot(222)
ax.plot(x, f1, '-', color='gray')
ax.fill_between(x, f1 - 2 * f1_err, f1 + 2 * f1_err, color='gray', alpha=0.3)
ax.plot(x1, y1, '.k', ms=6)
#ax.plot(x,tmp,'r')
lines2 = ax.plot(x, f1draws.T, lw=4, alpha=0.03+0.3/Ndraws, color='blue')
plugins.connect(fig, HighlightLines(lines2))
# third: plot a constrained function with errors
ax = fig.add_subplot(223)
ax.plot(x, f2, '-', color='gray')
ax.fill_between(x, f2 - 2 * f2_err, f2 + 2 * f2_err, color='gray', alpha=0.3)
ax.errorbar(x1, y1, dy2, fmt='.k', ms=6)
lines3 = ax.plot(x, f2draws.T, lw=4, alpha=0.03+0.3/Ndraws, color='blue')
plugins.connect(fig, HighlightLines(lines3))
ax.set_xlabel('$x$')
ax.set_ylabel('$f(x)$')
# third: plot a more constrained function with errors
ax = fig.add_subplot(224)
ax.plot(x, f3, '-', color='gray')
ax.fill_between(x, f3 - 2 * f3_err, f3 + 2 * f3_err, color='gray', alpha=0.3)
ax.errorbar(x3, y3, dy3, fmt='.k', ms=6)
ax.plot(x, np.cos(x), ':k')
ax.set_xlabel('$x$')
for ax in fig.axes:
ax.set_xlim(0, 10)
#plt.show()
mpld3.display()
Out[12]:
In [7]:
In [7]: