.. Copyright (C) 2014 Jaakko Luttinen This file is licensed under the MIT License. See LICENSE for a text of the license. .. testsetup:: # This is the PCA model from the previous sections import numpy as np np.random.seed(1) from bayespy.nodes import GaussianARD, Gamma, Dot D = 3 X = GaussianARD(0, 1, shape=(D,), plates=(1,100), name='X') alpha = Gamma(1e-3, 1e-3, plates=(D,), name='alpha') C = GaussianARD(0, alpha, shape=(D,), plates=(10,1), name='C') F = Dot(C, X) tau = Gamma(1e-3, 1e-3, name='tau') Y = GaussianARD(F, tau, name='Y') c = np.random.randn(10, 2) x = np.random.randn(2, 100) data = np.dot(c, x) + 0.1*np.random.randn(10, 100) Y.observe(data) Y.observe(data, mask=[[True], [False], [False], [True], [True], [False], [True], [True], [True], [False]]) X.initialize_from_parameters(np.random.randn(1, 100, D), 10) from bayespy.inference import VB Q = VB(Y, C, X, alpha, tau) X.initialize_from_parameters(np.random.randn(1, 100, D), 10) from bayespy.inference.vmp import transformations rotX = transformations.RotateGaussianARD(X) rotC = transformations.RotateGaussianARD(C, alpha) R = transformations.RotationOptimizer(rotC, rotX, D) Q = VB(Y, C, X, alpha, tau) Q.callback = R.rotate Q.update(repeat=1000, tol=1e-6, verbose=False) from bayespy.nodes import Gaussian Examining the results ===================== After the results have been obtained, it is important to be able to examine the results easily. The results can be examined either numerically by inspecting numerical arrays or visually by plotting distributions of the nodes. In addition, the posterior distributions can be visualized during the learning algorithm and the results can saved into a file. Plotting the results -------------------- .. currentmodule:: bayespy The module :mod:`plot` offers some plotting basic functionality: >>> import bayespy.plot as bpplt .. currentmodule:: bayespy.plot The module contains ``matplotlib.pyplot`` module if the user needs that. For instance, interactive plotting can be enabled as: >>> bpplt.pyplot.ion() The :mod:`plot` module contains some functions but it is not a very comprehensive collection, thus the user may need to write some problem- or model-specific plotting functions. The current collection is: * :func:`pdf`: show probability density function of a scalar * :func:`contour`: show probability density function of two-element vector * :func:`hinton`: show the Hinton diagram * :func:`plot`: show value as a function The probability density function of a scalar random variable can be plotted using the function :func:`pdf`: >>> bpplt.pyplot.figure() >>> bpplt.pdf(Q['tau'], np.linspace(60, 140, num=100)) [] .. plot:: # This is the PCA model from the previous sections import numpy as np np.random.seed(1) from bayespy.nodes import GaussianARD, Gamma, Dot D = 3 X = GaussianARD(0, 1, shape=(D,), plates=(1,100), name='X') alpha = Gamma(1e-3, 1e-3, plates=(D,), name='alpha') C = GaussianARD(0, alpha, shape=(D,), plates=(10,1), name='C') F = Dot(C, X) tau = Gamma(1e-3, 1e-3, name='tau') Y = GaussianARD(F, tau) c = np.random.randn(10, 2) x = np.random.randn(2, 100) data = np.dot(c, x) + 0.1*np.random.randn(10, 100) Y.observe(data) Y.observe(data, mask=[[True], [False], [False], [True], [True], [False], [True], [True], [True], [False]]) from bayespy.inference import VB Q = VB(Y, C, X, alpha, tau) X.initialize_from_parameters(np.random.randn(1, 100, D), 10) from bayespy.inference.vmp import transformations rotX = transformations.RotateGaussianARD(X) rotC = transformations.RotateGaussianARD(C, alpha) R = transformations.RotationOptimizer(rotC, rotX, D) Q = VB(Y, C, X, alpha, tau) Q.callback = R.rotate Q.update(repeat=1000, tol=1e-6) Q.update(repeat=50, tol=np.nan) import bayespy.plot as bpplt bpplt.pdf(Q['tau'], np.linspace(60, 140, num=100)) The variable ``tau`` models the inverse variance of the noise, for which the true value is :math:`0.1^{-2}=100`. Thus, the posterior captures the true value quite accurately. Similarly, the function :func:`contour` can be used to plot the probability density function of a 2-dimensional variable, for instance: >>> V = Gaussian([3, 5], [[4, 2], [2, 5]]) >>> bpplt.pyplot.figure() >>> bpplt.contour(V, np.linspace(1, 5, num=100), np.linspace(3, 7, num=100)) .. plot:: # This is the PCA model from the previous sections import numpy as np np.random.seed(1) from bayespy.nodes import GaussianARD, Gamma, Dot D = 3 X = GaussianARD(0, 1, shape=(D,), plates=(1,100), name='X') alpha = Gamma(1e-3, 1e-3, plates=(D,), name='alpha') C = GaussianARD(0, alpha, shape=(D,), plates=(10,1), name='C') F = Dot(C, X) tau = Gamma(1e-3, 1e-3, name='tau') Y = GaussianARD(F, tau) c = np.random.randn(10, 2) x = np.random.randn(2, 100) data = np.dot(c, x) + 0.1*np.random.randn(10, 100) Y.observe(data) Y.observe(data, mask=[[True], [False], [False], [True], [True], [False], [True], [True], [True], [False]]) from bayespy.inference import VB Q = VB(Y, C, X, alpha, tau) X.initialize_from_parameters(np.random.randn(1, 100, D), 10) from bayespy.inference.vmp import transformations rotX = transformations.RotateGaussianARD(X) rotC = transformations.RotateGaussianARD(C, alpha) R = transformations.RotationOptimizer(rotC, rotX, D) Q = VB(Y, C, X, alpha, tau) Q.callback = R.rotate Q.update(repeat=1000, tol=1e-6) Q.update(repeat=50, tol=np.nan) import bayespy.plot as bpplt from bayespy.nodes import Gaussian V = Gaussian([3, 5], [[4, 2], [2, 5]]) bpplt.contour(V, np.linspace(1, 5, num=100), np.linspace(3, 7, num=100)) Both :func:`pdf` and :func:`contour` require that the user provides the grid on which the probability density function is computed. They also support several keyword arguments for modifying the output, similarly as ``plot`` and ``contour`` in ``matplotlib.pyplot``. These functions can be used only for stochastic nodes. A few other plot types are also available as built-in functions. A Hinton diagram can be plotted as: >>> bpplt.pyplot.figure() >>> bpplt.hinton(C) .. plot:: # This is the PCA model from the previous sections import numpy as np np.random.seed(1) from bayespy.nodes import GaussianARD, Gamma, Dot D = 3 X = GaussianARD(0, 1, shape=(D,), plates=(1,100), name='X') alpha = Gamma(1e-3, 1e-3, plates=(D,), name='alpha') C = GaussianARD(0, alpha, shape=(D,), plates=(10,1), name='C') F = Dot(C, X) tau = Gamma(1e-3, 1e-3, name='tau') Y = GaussianARD(F, tau) c = np.random.randn(10, 2) x = np.random.randn(2, 100) data = np.dot(c, x) + 0.1*np.random.randn(10, 100) Y.observe(data) Y.observe(data, mask=[[True], [False], [False], [True], [True], [False], [True], [True], [True], [False]]) from bayespy.inference import VB Q = VB(Y, C, X, alpha, tau) X.initialize_from_parameters(np.random.randn(1, 100, D), 10) from bayespy.inference.vmp import transformations rotX = transformations.RotateGaussianARD(X) rotC = transformations.RotateGaussianARD(C, alpha) R = transformations.RotationOptimizer(rotC, rotX, D) Q = VB(Y, C, X, alpha, tau) Q.callback = R.rotate Q.update(repeat=1000, tol=1e-6) Q.update(repeat=50, tol=np.nan) import bayespy.plot as bpplt from bayespy.nodes import Gaussian bpplt.hinton(C) The diagram shows the elements of the matrix :math:`C`. The size of the filled rectangle corresponds to the absolute value of the element mean, and white and black correspond to positive and negative values, respectively. The non-filled rectangle shows standard deviation. From this diagram it is clear that the third column of :math:`C` has been pruned out and the rows that were missing in the data have zero mean and column-specific variance. The function :func:`hinton` is a simple wrapper for node-specific Hinton diagram plotters, such as :func:`gaussian_hinton` and :func:`dirichlet_hinton`. Thus, the keyword arguments depend on the node which is plotted. Another plotting function is :func:`plot`, which just plots the values of the node over one axis as a function: >>> bpplt.pyplot.figure() >>> bpplt.plot(X, axis=-2) .. plot:: # This is the PCA model from the previous sections import numpy as np np.random.seed(1) from bayespy.nodes import GaussianARD, Gamma, Dot D = 3 X = GaussianARD(0, 1, shape=(D,), plates=(1,100), name='X') alpha = Gamma(1e-3, 1e-3, plates=(D,), name='alpha') C = GaussianARD(0, alpha, shape=(D,), plates=(10,1), name='C') F = Dot(C, X) tau = Gamma(1e-3, 1e-3, name='tau') Y = GaussianARD(F, tau) c = np.random.randn(10, 2) x = np.random.randn(2, 100) data = np.dot(c, x) + 0.1*np.random.randn(10, 100) Y.observe(data) Y.observe(data, mask=[[True], [False], [False], [True], [True], [False], [True], [True], [True], [False]]) from bayespy.inference import VB Q = VB(Y, C, X, alpha, tau) X.initialize_from_parameters(np.random.randn(1, 100, D), 10) from bayespy.inference.vmp import transformations rotX = transformations.RotateGaussianARD(X) rotC = transformations.RotateGaussianARD(C, alpha) R = transformations.RotationOptimizer(rotC, rotX, D) Q = VB(Y, C, X, alpha, tau) Q.callback = R.rotate Q.update(repeat=1000, tol=1e-6) Q.update(repeat=50, tol=np.nan) import bayespy.plot as bpplt from bayespy.nodes import Gaussian bpplt.plot(X, axis=-2) Now, the ``axis`` is the second last axis which corresponds to :math:`n=0,\ldots,N-1`. As :math:`D=3`, there are three subplots. For Gaussian variables, the function shows the mean and two standard deviations. The plot shows that the third component has been pruned out, thus the method has been able to recover the true dimensionality of the latent space. It also has similar keyword arguments to ``plot`` function in ``matplotlib.pyplot``. Again, :func:`plot` is a simple wrapper over node-specific plotting functions, thus it supports only some node classes. Monitoring during the inference algorithm ----------------------------------------- It is possible to plot the distribution of the nodes during the learning algorithm. This is useful when the user is interested to see how the distributions evolve during learning and what is happening to the distributions. In order to utilize monitoring, the user must set plotters for the nodes that he or she wishes to monitor. This can be done either when creating the node or later at any time. The plotters are set by creating a plotter object and providing this object to the node. The plotter is a wrapper of one of the plotting functions mentioned above: :class:`PDFPlotter`, :class:`ContourPlotter`, :class:`HintonPlotter` or :class:`FunctionPlotter`. Thus, our example model could use the following plotters: >>> tau.set_plotter(bpplt.PDFPlotter(np.linspace(60, 140, num=100))) >>> C.set_plotter(bpplt.HintonPlotter()) >>> X.set_plotter(bpplt.FunctionPlotter(axis=-2)) These could have been given at node creation as a keyword argument ``plotter``: >>> V = Gaussian([3, 5], [[4, 2], [2, 5]], ... plotter=bpplt.ContourPlotter(np.linspace(1, 5, num=100), ... np.linspace(3, 7, num=100))) When the plotter is set, one can use the ``plot`` method of the node to perform plotting: >>> V.plot() Nodes can also be plotted using the ``plot`` method of the inference engine: >>> Q.plot('C') This method remembers the figure in which a node has been plotted and uses that every time it plots the same node. In order to monitor the nodes during learning, it is possible to use the keyword argument ``plot``: >>> Q.update(repeat=5, plot=True, tol=np.nan) Iteration 19: loglike=-1.221354e+02 (... seconds) Iteration 20: loglike=-1.221354e+02 (... seconds) Iteration 21: loglike=-1.221354e+02 (... seconds) Iteration 22: loglike=-1.221354e+02 (... seconds) Iteration 23: loglike=-1.221354e+02 (... seconds) Each node which has a plotter set will be plotted after it is updated. Note that this may slow down the inference significantly if the plotting operation is time consuming. Posterior parameters and moments -------------------------------- If the built-in plotting functions are not sufficient, it is possible to use ``matplotlib.pyplot`` for custom plotting. Each node has ``get_moments`` method which returns the moments and they can be used for plotting. Stochastic exponential family nodes have natural parameter vectors which can also be used. In addition to plotting, it is also possible to just print the moments or parameters in the console. Saving and loading results -------------------------- .. currentmodule:: bayespy.inference The results of the inference engine can be easily saved and loaded using :func:`VB.save` and :func:`VB.load` methods: >>> import tempfile >>> filename = tempfile.mkstemp(suffix='.hdf5')[1] >>> Q.save(filename=filename) >>> Q.load(filename=filename) The results are stored in a HDF5 file. The user may set an autosave file in which the results are automatically saved regularly. Autosave filename can be set at creation time by ``autosave_filename`` keyword argument or later using :func:`VB.set_autosave` method. If autosave file has been set, the :func:`VB.save` and :func:`VB.load` methods use that file by default. In order for the saving to work, all stochastic nodes must have been given (unique) names. However, note that these methods do *not* save nor load the node definitions. It means that the user must create the nodes and the inference engine and then use :func:`VB.load` to set the state of the nodes and the inference engine. If there are any differences in the model that was saved and the one which is tried to update using loading, then loading does not work. Thus, the user should keep the model construction unmodified in a Python file in order to be able to load the results later. Or if the user wishes to share the results, he or she must share the model construction Python file with the HDF5 results file.