################################################################################
# Copyright (C) 2011-2013 Jaakko Luttinen
#
# This file is licensed under the MIT License.
################################################################################
"""
Functions for plotting nodes.
Functions
=========
.. currentmodule:: bayespy.plot
.. autosummary::
:toctree: generated/
pdf
contour
plot
hinton
gaussian_mixture_2d
ellipse_from_cov
ellipse_from_precision
Plotters
========
.. autosummary::
:toctree: generated/
Plotter
PDFPlotter
ContourPlotter
HintonPlotter
FunctionPlotter
GaussianTimeseriesPlotter
CategoricalMarkovChainPlotter
"""
import os, sys
############################################################################
# A STUPID WORKAROUND FOR A MATPLOTLIB 1.4.0 BUG RELATED TO INTERACTIVE MODE
# See: https://github.com/matplotlib/matplotlib/issues/3505
import __main__
if hasattr(__main__, '__file__'):
sys.ps1 = ('WORKAROUND FOR A BUG #3505 IN MATPLOTLIB.\n'
'IF YOU SEE THIS MESSAGE, TRY MATPLOTLIB!=1.4.0.')
# This workaround does not work on Python shell, only on stand-alone scripts
# and IPython. A better solution: require MPL!=1.4.0.
#############################################################################
import numpy as np
import scipy.sparse as sp
import scipy
from scipy import special
import matplotlib.pyplot as plt
from matplotlib import animation
from matplotlib import colors
from matplotlib import patches
#from matplotlib.pyplot import *
from bayespy.inference.vmp.nodes.categorical import CategoricalMoments
from bayespy.inference.vmp.nodes.gaussian import (GaussianMoments,
GaussianWishartMoments)
from bayespy.inference.vmp.nodes.beta import BetaMoments
from bayespy.inference.vmp.nodes.beta import DirichletMoments
from bayespy.inference.vmp.nodes.bernoulli import BernoulliMoments
from bayespy.inference.vmp.nodes.categorical import CategoricalMoments
from bayespy.inference.vmp.nodes.gamma import GammaMoments
from bayespy.inference.vmp.nodes.node import Node, Moments
from bayespy.utils import (misc,
random,
linalg)
# Users can use pyplot via this module
import matplotlib
mpl = matplotlib
pyplot = plt
def interactive(function):
"""A decorator for forcing functions to use the interactive mode.
Parameters
----------
function : callable
The function to be decorated
"""
def new_function(*args, **kwargs):
if mpl.is_interactive():
was_interactive = True
else:
was_interactive = False
mpl.interactive(True)
retval = function(*args, **kwargs)
if not was_interactive:
mpl.interactive(False)
return retval
return new_function
def _subplots(plotfunc, *args, fig=None, kwargs=None):
"""Create a collection of subplots
Each subplot is created with the same plotting function.
Inputs are given as pairs:
(x, 3), (y, 2), ...
where x,y,... are the input arrays and 3,2,... are the ndim
parameters. The last ndim axes of each array are interpreted as a
single element to the plotting function.
All high-level plotting functions should wrap low-level plotting
functions with this function in order to generate subplots for
plates.
"""
if kwargs is None:
kwargs = {}
if fig is None:
fig = plt.gcf()
# Parse shape and plates of each input array
shapes = [np.shape(x)[-n:] if n > 0 else ()
for (x,n) in args]
plates = [np.shape(x)[:-n] if n > 0 else np.shape(x)
for (x,n) in args]
# Get the full grid shape of the subplots
broadcasted_plates = misc.broadcasted_shape(*plates)
# Subplot indexing layout
M = int(np.prod(broadcasted_plates[-2::-2]))
N = int(np.prod(broadcasted_plates[-1::-2]))
strides_subplot = [np.prod(broadcasted_plates[(j+2)::2]) * N
if ((len(broadcasted_plates)-j) % 2) == 0 else
np.prod(broadcasted_plates[(j+2)::2])
for j in range(len(broadcasted_plates))]
# Plot each subplot
for ind in misc.nested_iterator(broadcasted_plates):
# Get the list of inputs for this subplot
broadcasted_args = []
for n in range(len(args)):
i = misc.safe_indices(ind, plates[n])
broadcasted_args.append(args[n][0][i])
# Plot the subplot using the given function
ind_subplot = int(np.einsum('i,i', ind, strides_subplot))
axes = fig.add_subplot(M, N, ind_subplot+1)
plotfunc(*broadcasted_args, axes=axes, **kwargs)
[docs]def pdf(Z, x, *args, name=None, axes=None, fig=None, **kwargs):
"""
Plot probability density function of a scalar variable.
Parameters
----------
Z : node or function
Stochastic node or log pdf function
x : array
Grid points
"""
# TODO: Make it possible to plot a plated variable using _subplots function.
if axes is None and fig is None:
axes = plt.gca()
else:
if fig is None:
fig = plt.gcf()
axes = fig.add_subplot(111)
try:
lpdf = Z.logpdf(x)
except AttributeError:
lpdf = Z(x)
p = np.exp(lpdf)
retval = axes.plot(x, p, *args, **kwargs)
if name is None:
try:
name = Z.name
except AttributeError:
pass
if name:
axes.set_title(r'$q(%s)$' % (name))
axes.set_xlabel(r'$%s$' % (name))
return retval
[docs]def contour(Z, x, y, n=None, axes=None, fig=None, **kwargs):
"""
Plot 2-D probability density function of a 2-D variable.
Parameters
----------
Z : node or function
Stochastic node or log pdf function
x : array
Grid points on x axis
y : array
Grid points on y axis
"""
# TODO: Make it possible to plot a plated variable using _subplots function.
if axes is None and fig is None:
axes = plt.gca()
else:
if fig is None:
fig = plt.gcf()
axes = fig.add_subplot(111)
XY = misc.grid(x, y)
try:
lpdf = Z.logpdf(XY)
except AttributeError:
lpdf = Z(XY)
p = np.exp(lpdf)
shape = (np.size(x), np.size(y))
X = np.reshape(XY[:,0], shape)
Y = np.reshape(XY[:,1], shape)
P = np.reshape(p, shape)
if n is not None:
levels = np.linspace(0, np.amax(P), num=n+2)[1:-1]
return axes.contour(X, Y, P, levels, **kwargs)
else:
return axes.contour(X, Y, P, **kwargs)
def plot_gaussian_mc(X, scale=2, **kwargs):
"""
Plot Gaussian Markov chain as a 1-D function
Parameters
----------
X : node
Node with Gaussian Markov chain moments.
"""
timeseries_gaussian(X, axis=-2, scale=scale, **kwargs)
def plot_bernoulli(X, axis=-1, scale=2, **kwargs):
"""
Plot Bernoulli node as a 1-D function
"""
X = X._ensure_moments(X, BernoulliMoments)
u_X = X.get_moments()
z = u_X[0]
return _timeseries_mean_and_error(z, None, axis=axis, **kwargs)
def plot_gaussian(X, axis=-1, scale=2, **kwargs):
"""
Plot Gaussian node as a 1-D function
Parameters
----------
X : node
Node with Gaussian moments.
axis : int
The index of the time axis.
"""
X = X._ensure_moments(X, GaussianMoments, ndim=0)
u_X = X.get_moments()
x = u_X[0]
xx = misc.get_diag(u_X[1], ndim=len(X.dims[0]))
std = scale * np.sqrt(xx - x**2)
#std = scale * np.sqrt(np.einsum('...ii->...i', xx) - x**2)
return _timeseries_mean_and_error(x, std, axis=axis, **kwargs)
[docs]def plot(Y, axis=-1, scale=2, center=False, **kwargs):
"""
Plot a variable or an array as 1-D function with errorbars
"""
if misc.is_numeric(Y):
return _timeseries_mean_and_error(Y, None, axis=axis, center=center, **kwargs)
if isinstance(Y, Node):
# Try Bernoulli plotting
try:
Y = Y._ensure_moments(Y, BernoulliMoments)
except BernoulliMoments.NoConverterError:
pass
else:
return plot_bernoulli(Y, axis=axis, scale=scale, center=center, **kwargs)
# Try Gaussian plotting
try:
Y = Y._ensure_moments(Y, GaussianMoments, ndim=0)
except GaussianMoments.NoConverterError:
pass
else:
return plot_gaussian(Y, axis=axis, scale=scale, center=center, **kwargs)
(mu, var) = Y.get_mean_and_variance()
std = np.sqrt(var)
return _timeseries_mean_and_error(mu, std,
axis=axis,
scale=scale,
center=center,
**kwargs)
# Some backward compatibility
def timeseries_gaussian_mc(*args, center=True, **kwargs):
return plot_gaussian_mc(*args, center=center, **kwargs)
def timeseries_gaussian(*args, center=True, **kwargs):
return plot_gaussian(*args, center=center, **kwargs)
timeseries_normal = timeseries_gaussian
def timeseries(*args, center=True, **kwargs):
return plot(*args, center=center, **kwargs)
def _timeseries_mean_and_error(y, std, *args, axis=-1, center=True, fig=None, axes=None, **kwargs):
# TODO/FIXME: You must multiply by ones(plates) in order to plot
# broadcasted plates properly
if fig is None:
fig = plt.gcf()
y = np.atleast_1d(y)
shape = list(np.shape(y))
# Get and remove the length of the time axis
T = shape.pop(axis)
# Move time axis to first
y = np.rollaxis(y, axis)
if std is not None:
std = np.rollaxis(std, axis)
y = np.reshape(y, (T, -1))
if std is not None:
std = np.reshape(std, (T, -1))
# Remove 1s
shape = [s for s in shape if s > 1]
# Calculate number of rows and columns
shape = misc.multiply_shapes(shape, (1,1))
if len(shape) > 2:
raise Exception("Can plot only in 2 dimensions (rows and columns)")
(M, N) = shape
# Prefer plotting to rows
if M == 1:
M = N
N = 1
# Plot each timeseries
if axes is None:
ax0 = fig.add_subplot(M, N, 1)
for i in range(M*N):
if axes is None:
if i > 0:
# Share x axis between all subplots
ax = fig.add_subplot(M, N, i+1, sharex=ax0)
else:
ax = ax0
# Autoscale the axes to data and use tight y and x axes
ax.autoscale(enable=True, tight=True)
ax.set_ylim(auto=True)
if i < (M-1)*N:
# Remove x tick labels from other than the last row
plt.setp(ax.get_xticklabels(), visible=False)
else:
ax = axes[i]
if std is None:
errorplot(y=y[:,i], axes=ax, **kwargs)
else:
if len(args) > 0:
raise Exception("Can't handle extra arguments")
errorplot(y=y[:,i], error=std[:,i], axes=ax, **kwargs)
if center:
# Center the zero level on y-axis
ylim = ax.get_ylim()
vmax = np.max(np.abs(ylim))
ax.set_ylim([-vmax, vmax])
if axes is None:
# Remove height space between subplots
fig.subplots_adjust(hspace=0)
def _blob(axes, x, y, area, colour):
"""
Draws a square-shaped blob with the given area (< 1) at
the given coordinates.
"""
hs = np.sqrt(area) / 2
xcorners = np.array([x - hs, x + hs, x + hs, x - hs])
ycorners = np.array([y - hs, y - hs, y + hs, y + hs])
axes.fill(xcorners, ycorners, colour, edgecolor=colour)
def _rectangle(axes, x, y, width, height, **kwargs):
_x = x - width/2
_y = y - height/2
rectangle = plt.Rectangle((_x, _y),
width,
height,
**kwargs)
axes.add_patch(rectangle)
return
[docs]def gaussian_mixture_2d(X, alpha=None, scale=2, fill=False, axes=None, **kwargs):
"""
Plot Gaussian mixture as ellipses in 2-D
Parameters
----------
X : Mixture node
alpha : Dirichlet-like node (optional)
Probabilities for the clusters
scale : float (optional)
Scale for the covariance ellipses (by default, 2)
"""
if axes is None:
axes = plt.gca()
mu_Lambda = X._ensure_moments(X.parents[1], GaussianWishartMoments)
(mu, _, Lambda, _) = mu_Lambda.get_moments()
mu = np.linalg.solve(Lambda, mu)
if len(mu_Lambda.plates) != 1:
raise NotImplementedError("Not yet implemented for more plates")
K = mu_Lambda.plates[0]
width = np.zeros(K)
height = np.zeros(K)
angle = np.zeros(K)
for k in range(K):
m = mu[k]
L = Lambda[k]
(u, W) = scipy.linalg.eigh(L)
u[0] = np.sqrt(1/u[0])
u[1] = np.sqrt(1/u[1])
width[k] = 2*u[0]
height[k] = 2*u[1]
angle[k] = np.arctan(W[0,1] / W[0,0])
angle = 180 * angle / np.pi
mode_height = 1 / (width * height)
# Use cluster probabilities to adjust alpha channel
if alpha is not None:
# Compute the normalized probabilities in a numerically stable way
logsum_p = misc.logsumexp(alpha.u[0], axis=-1, keepdims=True)
logp = alpha.u[0] - logsum_p
p = np.exp(logp)
# Visibility is based on cluster mode peak height
visibility = mode_height * p
visibility /= np.amax(visibility)
else:
visibility = np.ones(K)
for k in range(K):
ell = mpl.patches.Ellipse(mu[k], scale*width[k], scale*height[k],
angle=(180+angle[k]),
fill=fill,
alpha=visibility[k],
**kwargs)
axes.add_artist(ell)
plt.axis('equal')
# If observed, plot the data too
if np.any(X.observed):
mask = np.array(X.observed) * np.ones(X.plates, dtype=np.bool)
y = X.u[0][mask]
plt.plot(y[:,0], y[:,1], 'r.')
return
def _hinton(W, error=None, vmax=None, square=False, axes=None):
"""
Draws a Hinton diagram for visualizing a weight matrix.
Temporarily disables matplotlib interactive mode if it is on,
otherwise this takes forever.
Originally copied from
http://wiki.scipy.org/Cookbook/Matplotlib/HintonDiagrams
"""
if axes is None:
axes = plt.gca()
W = misc.atleast_nd(W, 2)
(height, width) = W.shape
if not vmax:
#vmax = 2**np.ceil(np.log(np.max(np.abs(W)))/np.log(2))
if error is not None:
vmax = np.max(np.abs(W) + error)
else:
vmax = np.max(np.abs(W))
axes.fill(0.5+np.array([0,width,width,0]),
0.5+np.array([0,0,height,height]),
'gray')
if square:
axes.set_aspect('equal')
axes.set_ylim(0.5, height+0.5)
axes.set_xlim(0.5, width+0.5)
axes.set_xticks([])
axes.set_yticks([])
axes.invert_yaxis()
for x in range(width):
for y in range(height):
_x = x+1
_y = y+1
w = W[y,x]
_w = np.abs(w)
if w > 0:
_c = 'white'
else:
_c = 'black'
if error is not None:
e = error[y,x]
if e < 0:
print(e, _w, vmax)
raise Exception("BUG? Negative error")
if _w + e > vmax:
print(e, _w, vmax)
raise Exception("BUG? Value+error greater than max")
_rectangle(axes,
_x,
_y,
min(1, np.sqrt((_w+e)/vmax)),
min(1, np.sqrt((_w+e)/vmax)),
edgecolor=_c,
fill=False)
_blob(axes, _x, _y, min(1, _w/vmax), _c)
def matrix(A, axes=None):
if axes is None:
axes = plt.gca()
A = np.atleast_2d(A)
vmax = np.max(np.abs(A))
return axes.imshow(A,
interpolation='nearest',
cmap='RdBu_r',
vmin=-vmax,
vmax=vmax)
def new_matrix(A, vmax=None):
A = np.atleast_2d(A)
if vmax is None:
vmax = np.max(np.abs(A))
(M, N) = np.shape(A)
for i in range(M):
for j in range(N):
pass
def gaussian_hinton(X, rows=None, cols=None, scale=1, fig=None):
"""
Plot the Hinton diagram of a Gaussian node
"""
if fig is None:
fig = plt.gcf()
# Get mean and second moment
X = X._ensure_moments(X, GaussianMoments, ndim=0)
(x, xx) = X.get_moments()
ndim = len(X.dims[0])
shape = X.get_shape(0)
size = len(X.get_shape(0))
# Compute standard deviation
xx = misc.get_diag(xx, ndim=ndim)
std = np.sqrt(xx - x**2)
# Force explicit elements when broadcasting
x = x * np.ones(shape)
std = std * np.ones(shape)
if rows is None:
rows = np.nan
if cols is None:
cols = np.nan
# Preprocess the axes to 0,...,ndim
if rows < 0:
rows += size
if cols < 0:
cols += size
if rows < 0 or rows >= size:
raise ValueError("Row axis invalid")
if cols < 0 or cols >= size:
raise ValueError("Column axis invalid")
# Remove non-row and non-column axes that have length 1
squeezed_shape = list(shape)
for i in reversed(range(len(shape))):
if shape[i] == 1 and i != rows and i != cols:
squeezed_shape.pop(i)
if i < cols:
cols -= 1
if i < rows:
rows -= 1
x = np.reshape(x, squeezed_shape)
std = np.reshape(std, squeezed_shape)
if np.ndim(x) < 2:
cols += 2 - np.ndim(x)
rows += 2 - np.ndim(x)
x = np.atleast_2d(x)
std = np.atleast_2d(std)
size = np.ndim(x)
if np.isnan(cols):
if rows != size - 1:
cols = size - 1
else:
cols = size - 2
if np.isnan(rows):
if cols != size - 1:
rows = size - 1
else:
rows = size - 2
# Put the row and column axes to the end
axes = [i for i in range(size) if i not in (rows, cols)] + [rows, cols]
x = np.transpose(x, axes=axes)
std = np.transpose(std, axes=axes)
vmax = np.max(np.abs(x) + scale*std)
if scale == 0:
_subplots(_hinton, (x, 2), fig=fig, kwargs=dict(vmax=vmax))
else:
def plotfunc(z, e, **kwargs):
return _hinton(z, error=e, **kwargs)
_subplots(plotfunc, (x, 2), (scale*std, 2), fig=fig, kwargs=dict(vmax=vmax))
def _hinton_figure(x, rows=None, cols=None, fig=None, square=True):
"""
Plot the Hinton diagram of a Gaussian node
"""
scale = 0
std = 0
if fig is None:
fig = plt.gcf()
# Get mean and second moment
shape = np.shape(x)
size = np.ndim(x)
if rows is None:
rows = np.nan
if cols is None:
cols = np.nan
# Preprocess the axes to 0,...,ndim
if rows < 0:
rows += size
if cols < 0:
cols += size
if rows < 0 or rows >= size:
raise ValueError("Row axis invalid")
if cols < 0 or cols >= size:
raise ValueError("Column axis invalid")
# Remove non-row and non-column axes that have length 1
squeezed_shape = list(shape)
for i in reversed(range(len(shape))):
if shape[i] == 1 and i != rows and i != cols:
squeezed_shape.pop(i)
if i < cols:
cols -= 1
if i < rows:
rows -= 1
x = np.reshape(x, squeezed_shape)
size = np.ndim(x)
if np.isnan(cols):
if rows != size - 1:
cols = size - 1
else:
cols = size - 2
if np.isnan(rows):
if cols != size - 1:
rows = size - 1
else:
rows = size - 2
# Put the row and column axes to the end
if np.ndim(x) >= 2:
axes = [i for i in range(size) if i not in (rows, cols)] + [rows, cols]
x = np.transpose(x, axes=axes)
#std = np.transpose(std, axes=axes)
vmax = np.max(np.abs(x) + scale*std)
kw = dict(vmax=vmax, square=square)
if scale == 0:
_subplots(_hinton, (x, 2), fig=fig, kwargs=kw)
else:
def plotfunc(z, e, **kwargs):
return _hinton(z, error=e, **kwargs)
_subplots(plotfunc, (x, 2), (scale*std, 2), fig=fig, kwargs=kw)
# For backwards compatibility:
gaussian_array = gaussian_hinton
def timeseries_categorical_mc(Z, fig=None):
if fig is None:
fig = plt.gcf()
# Make sure that the node is categorical
Z = Z._ensure_moments(Z, CategoricalMoments, categories=None)
# Get expectations (and broadcast explicitly)
z = Z._message_to_child()[0] * np.ones(Z.get_shape(0))
# Compute the subplot layout
z = misc.atleast_nd(z, 4)
if np.ndim(z) != 4:
raise ValueError("Can not plot arrays with over 4 axes")
M = np.shape(z)[0]
N = np.shape(z)[1]
# Plot Hintons
for i in range(M):
for j in range(N):
axes = fig.add_subplot(M, N, i*N+j+1)
_hinton(z[i,j].T, vmax=1.0, square=False, axes=axes)
def gamma_hinton(alpha, square=True, **kwargs):
"""
Plot a beta distributed random variable as a Hinton diagram
"""
# Make sure that the node is beta
alpha = alpha._ensure_moments(alpha, GammaMoments)
# Compute exp( <log p> )
x = alpha.get_moments()[0]
# Explicit broadcasting
x = x * np.ones(alpha.plates)
# Plot Hinton diagram
return _hinton_figure(x, square=square, **kwargs)
def beta_hinton(P, square=True):
"""
Plot a beta distributed random variable as a Hinton diagram
"""
# Make sure that the node is beta
P = P._ensure_moments(P, BetaMoments)
# Compute exp( <log p> )
p = np.exp(P._message_to_child()[0][...,0])
# Explicit broadcasting
p = p * np.ones(P.plates)
# Plot Hinton diagram
return _hinton(p, vmax=1.0, square=square)
def dirichlet_hinton(P, square=True):
"""
Plot a beta distributed random variable as a Hinton diagram
"""
# Make sure that the node is beta
P = P._ensure_moments(P, DirichletMoments)
# Compute exp( <log p> )
p = np.exp(P._message_to_child()[0])
# Explicit broadcasting
p = p * np.ones(P.plates+(1,))
# Plot Hinton diagram
return _hinton(p, vmax=1.0, square=square)
def bernoulli_hinton(Z, square=True):
"""
Plot a Bernoulli distributed random variable as a Hinton diagram
"""
# Make sure that the node is Bernoulli
Z = Z._ensure_moments(Z, BernoulliMoments)
# Get <Z>
z = Z._message_to_child()[0]
# Explicit broadcasting
z = z * np.ones(Z.plates)
# Plot Hinton diagram
return _hinton(z, vmax=1.0, square=square)
def categorical_hinton(Z, square=True):
"""
Plot a Bernoulli distributed random variable as a Hinton diagram
"""
# Make sure that the node is Bernoulli
Z = Z._ensure_moments(Z, CategoricalMoments, categories=None)
# Get <Z>
z = Z._message_to_child()[0]
# Explicit broadcasting
z = z * np.ones(Z.plates+(1,))
# Plot Hinton diagram
return _hinton(np.squeeze(z), vmax=1.0, square=square)
[docs]def hinton(X, **kwargs):
r"""
Plot the Hinton diagram of a node
The keyword arguments depend on the node type. For some node types, the
diagram also shows uncertainty with non-filled rectangles. Currently,
beta-like, Gaussian-like and Dirichlet-like nodes are supported.
Parameters
----------
X : node
"""
if hasattr(X, "_ensure_moments"):
try:
X = X._ensure_moments(X, GaussianMoments, ndim=0)
except Moments.NoConverterError:
pass
else:
return gaussian_hinton(X, **kwargs)
try:
X = X._ensure_moments(X, GammaMoments)
except Moments.NoConverterError:
pass
else:
return gamma_hinton(X, **kwargs)
try:
X = X._ensure_moments(X, BetaMoments)
except Moments.NoConverterError:
pass
else:
return beta_hinton(X, **kwargs)
try:
X = X._ensure_moments(X, DirichletMoments)
except Moments.NoConverterError:
pass
else:
return dirichlet_hinton(X, **kwargs)
try:
X = X._ensure_moments(X, BernoulliMoments)
except Moments.NoConverterError:
pass
else:
return bernoulli_hinton(X, **kwargs)
try:
X = X._ensure_moments(X, CategoricalMoments, categories=None)
except Moments.NoConverterError:
pass
else:
return categorical_hinton(X, **kwargs)
return _hinton_figure(X, **kwargs)
[docs]class Plotter():
r"""
Wrapper for plotting functions and base class for node plotters
The purpose of this class is to collect all the parameters needed by a
plotting function and provide a callable interface which needs only the node
as the input.
Plotter instances are callable objects that plot a given node using a
specified plotting function.
Parameters
----------
plotter : function
Plotting function to use
args : defined by the plotting function
Additional inputs needed by the plotting function
kwargs : defined by the plotting function
Additional keyword arguments supported by the plotting function
Examples
--------
First, create a gamma variable:
>>> import numpy as np
>>> from bayespy.nodes import Gamma
>>> x = Gamma(4, 5)
The probability density function can be plotted as:
>>> import bayespy.plot as bpplt
>>> bpplt.pdf(x, np.linspace(0.1, 10, num=100)) # doctest: +ELLIPSIS
[<matplotlib.lines.Line2D object at 0x...>]
However, this can be problematic when one needs to provide a
plotting function for the inference engine as the inference engine
gives only the node as input. Thus, we need to create a simple
plotter wrapper:
>>> p = bpplt.Plotter(bpplt.pdf, np.linspace(0.1, 10, num=100))
Now, this callable object ``p`` needs only the node as the input:
>>> p(x) # doctest: +ELLIPSIS
[<matplotlib.lines.Line2D object at 0x...>]
Thus, it can be given to the inference engine to use as a plotting function:
>>> x = Gamma(4, 5, plotter=p)
>>> x.plot() # doctest: +ELLIPSIS
[<matplotlib.lines.Line2D object at 0x...>]
"""
[docs] def __init__(self, plotter, *args, **kwargs):
self._args = args
self._kwargs = kwargs
self._plotter = plotter
def __call__(self, X, fig=None, **kwargs):
"""
Plot the node using the specified plotting function
Parameters
----------
X : node
The plotted node
"""
kwargs_all = self._kwargs.copy()
kwargs_all.update(kwargs)
return self._plotter(X, *self._args, fig=fig, **kwargs_all)
[docs]class PDFPlotter(Plotter):
r"""
Plotter of probability density function of a scalar node
Parameters
----------
x_grid : array
Numerical grid on which the density function is computed and
plotted
See also
--------
pdf
"""
[docs] def __init__(self, x_grid, **kwargs):
super().__init__(pdf, x_grid, **kwargs)
[docs]class ContourPlotter(Plotter):
r"""
Plotter of probability density function of a two-dimensional node
Parameters
----------
x1_grid : array
Grid for the first dimension
x2_grid : array
Grid for the second dimension
See also
--------
contour
"""
[docs] def __init__(self, x1_grid, x2_grid, **kwargs):
super().__init__(contour, x1_grid, x2_grid, **kwargs)
[docs]class HintonPlotter(Plotter):
r"""
Plotter of the Hinton diagram of a node
See also
--------
hinton
"""
[docs] def __init__(self, **kwargs):
super().__init__(hinton, **kwargs)
[docs]class FunctionPlotter(Plotter):
r"""
Plotter of a node as a 1-dimensional function
See also
--------
plot
"""
[docs] def __init__(self, **kwargs):
super().__init__(plot, **kwargs)
class GaussianMarkovChainPlotter(Plotter):
r"""
Plotter of a Gaussian Markov chain as a timeseries
"""
def __init__(self, **kwargs):
super().__init__(timeseries_gaussian_mc, **kwargs)
[docs]class GaussianTimeseriesPlotter(Plotter):
r"""
Plotter of a Gaussian node as a timeseries
"""
[docs] def __init__(self, **kwargs):
super().__init__(timeseries_gaussian, **kwargs)
class GaussianHintonPlotter(Plotter):
r"""
Plotter of a Gaussian node as a Hinton diagram
"""
def __init__(self, **kwargs):
super().__init__(gaussian_array, **kwargs)
[docs]class CategoricalMarkovChainPlotter(Plotter):
r"""
Plotter of a Categorical timeseries
"""
[docs] def __init__(self, **kwargs):
super().__init__(timeseries_categorical_mc, **kwargs)
def matrix_animation(A, filename=None, fps=25, fig=None, **kwargs):
if fig is None:
fig = plt.gcf()
axes = fig.add_subplot(111)
A = np.atleast_3d(A)
vmax = np.max(np.abs(A))
x = axes.imshow(A[0],
interpolation='nearest',
cmap='RdBu_r',
vmin=-vmax,
vmax=vmax,
**kwargs)
s = axes.set_title('t = %d' % 0)
def animate(nframe):
s.set_text('t = %d' % nframe)
x.set_array(A[nframe])
return (x, s)
anim = animation.FuncAnimation(fig, animate,
frames=np.shape(A)[0],
interval=1000/fps,
blit=False,
repeat=False)
return anim
def save_animation(anim, filename, fps=25, bitrate=5000, fig=None):
# A bug in numpy/matplotlib causes this not to work in python3.3:
# https://github.com/matplotlib/matplotlib/issues/1891
#
# So the following command does not work currently..
#
# anim.save(filename, fps=fps)
if fig is None:
fig = plt.gcf()
writer = animation.FFMpegFileWriter(fps=fps, bitrate=bitrate)
writer.setup(fig, filename, 100)
anim.save(filename,
fps=fps,
writer=writer,
bitrate=bitrate)
return
def binary_matrix(A, axes=None):
if axes is None:
axes = plt.gca()
A = np.atleast_2d(A)
G = np.zeros(np.shape(A) + (3,))
G[A] = [0,0,0]
G[np.logical_not(A)] = [1,1,1]
axes.imshow(G, interpolation='nearest')
def gaussian_mixture_logpdf(x, w, mu, Sigma):
# Shape(x) = (N, D)
# Shape(w) = (K,)
# Shape(mu) = (K, D)
# Shape(Sigma) = (K, D, D)
# Shape(result) = (N,)
# Dimensionality
D = np.shape(x)[-1]
# Cholesky decomposition of the covariance matrix
U = linalg.chol(Sigma)
# Reshape x:
# Shape(x) = (N, 1, D)
x = np.expand_dims(x, axis=-2)
# (x-mu) and (x-mu)'*inv(Sigma)*(x-mu):
# Shape(v) = (N, K, D)
# Shape(z) = (N, K)
v = x - mu
z = np.einsum('...i,...i', v, linalg.chol_solve(U, v))
# Log-determinant of Sigma:
# Shape(ldet) = (K,)
ldet = linalg.chol_logdet(U)
# Compute log pdf for each cluster:
# Shape(lpdf) = (N, K)
lpdf = misc.gaussian_logpdf(z, 0, 0, ldet, D)
def matrixplot(A, colorbar=False, axes=None):
if axes is None:
axes = plt.gca()
if sp.issparse(A):
A = A.toarray()
axes.imshow(A, interpolation='nearest')
if colorbar:
plt.colorbar(ax=axes)
def contourplot(x1, x2, y, colorbar=False, filled=True, axes=None):
""" Plots 2D contour plot. x1 and x2 are 1D vectors, y contains
the function values. y.size must be x1.size*x2.size. """
if axes is None:
axes = plt.gca()
y = np.reshape(y, (len(x2),len(x1)))
if filled:
axes.contourf(x1, x2, y)
else:
axes.contour(x1, x2, y)
if colorbar:
plt.colorbar(ax=axes)
def errorplot(y=None, error=None, x=None, lower=None, upper=None,
color=(0,0,0,1), fillcolor=None, axes=None, **kwargs):
if axes is None:
axes = plt.gca()
# Default inputs
if x is None:
x = np.arange(np.size(y))
# Parse errors (lower=lower/error/upper, upper=upper/error/lower)
if lower is None:
if error is not None:
lower = error
elif upper is not None:
lower = upper
if upper is None:
if error is not None:
upper = error
elif lower is not None:
upper = lower
# Plot errors
if (lower is not None) and (upper is not None):
l = y - lower
u = y + upper
if fillcolor is None:
color = colors.ColorConverter().to_rgba(color)
fillcolor = tuple(color[:3]) + (0.2 * color[3],)
axes.fill_between(x,
l,
u,
facecolor=fillcolor,
edgecolor=(0, 0, 0, 0),
linewidth=1,
interpolate=True)
# Plot function
axes.plot(x, y, color=color, **kwargs)
def plotmatrix(X):
"""
Creates a matrix of marginal plots.
On diagonal, are marginal plots of each variable. Off-diagonal plot (i,j)
shows the joint marginal density of x_i and x_j.
"""
return X.plotmatrix()
def _pdf_t(mu, s2, nu, axes=None, scale=4, color='k'):
"""
"""
if axes is None:
axes = plt.gca()
s = np.sqrt(s2)
x = np.linspace(mu-scale*s, mu+scale*s, num=100)
y2 = (x-mu)**2 / s2
lpdf = random.t_logpdf(y2, np.log(s2), nu, 1)
p = np.exp(lpdf)
return axes.plot(x, p, color=color)
def _pdf_gamma(a, b, axes=None, scale=4, color='k'):
"""
"""
if axes is None:
axes = plt.gca()
if np.size(a) != 1 or np.size(b) != 1:
raise ValueError("Parameters must be scalars")
mean = a/b
v = scale*np.sqrt(a/b**2)
m = max(0, mean-v)
n = mean + v
x = np.linspace(m, n, num=100)
logx = np.log(x)
lpdf = random.gamma_logpdf(b*x,
logx,
a*logx,
a*np.log(b),
special.gammaln(a))
p = np.exp(lpdf)
return axes.plot(x, p, color=color)
def _contour_t(mu, Cov, nu, axes=None, scale=4, transpose=False, colors='k'):
"""
"""
if axes is None:
axes = plt.gca()
if np.shape(mu) != (2,) or np.shape(Cov) != (2,2) or np.shape(nu) != ():
print(np.shape(mu), np.shape(Cov), np.shape(nu))
raise ValueError("Only 2-d t-distribution allowed")
if transpose:
mu = mu[[1,0]]
Cov = Cov[np.ix_([1,0],[1,0])]
s = np.sqrt(np.diag(Cov))
x0 = np.linspace(mu[0]-scale*s[0], mu[0]+scale*s[0], num=100)
x1 = np.linspace(mu[1]-scale*s[1], mu[1]+scale*s[1], num=100)
X0X1 = misc.grid(x0, x1)
Y = X0X1 - mu
L = linalg.chol(Cov)
logdet_Cov = linalg.chol_logdet(L)
Z = linalg.chol_solve(L, Y)
Z = linalg.inner(Y, Z, ndim=1)
lpdf = random.t_logpdf(Z, logdet_Cov, nu, 2)
p = np.exp(lpdf)
shape = (np.size(x0), np.size(x1))
X0 = np.reshape(X0X1[:,0], shape)
X1 = np.reshape(X0X1[:,1], shape)
P = np.reshape(p, shape)
return axes.contour(X0, X1, P, colors=colors)
def _contour_gaussian_gamma(mu, s2, a, b, axes=None, transpose=False):
"""
"""
pass
[docs]def ellipse_from_cov(xy, cov, scale=2, **kwargs):
"""
Create an ellipse from a covariance matrix.
Parameters
----------
xy : np.ndarray
position of the ellipse
cov : np.ndarray
covariance matrix
scale : float
scale of the ellipse (default is three standard deviations)
kwargs : dict
keyword arguments passed on to `matplotlib.patches.Ellipse`
Returns
-------
ellipse : matplotlib.patches.Ellipse
"""
evals, evecs = np.linalg.eigh(cov)
angle = np.arctan2(*evecs[::-1, 0])
width, height = scale * np.sqrt(evals)
return patches.Ellipse(xy, width, height, np.rad2deg(angle), **kwargs)
[docs]def ellipse_from_precision(xy, precision, scale=2, **kwargs):
"""
Create an ellipse from a covariance matrix.
Parameters
----------
xy : np.ndarray
position of the ellipse
cov : np.ndarray
covariance matrix
scale : float
scale of the ellipse (default is three standard deviations)
kwargs : dict
keyword arguments passed on to `matplotlib.patches.Ellipse`
Returns
-------
ellipse : matplotlib.patches.Ellipse
"""
return ellipse_from_cov(xy, np.linalg.inv(precision), scale, **kwargs)