ENH: Add classes which define the sufficient statistics

These classes define the messages to children. Separating this protocol
from the nodes gains a few benefits:
  * Write statistic-specific features in one place only. For instance,
    covariance from Gaussian message.
  * Different nodes may have identically defined statistic so you need
    to implement related features only once. For instance, Gaussian and
    GaussianARD differ on the prior but the moments are the same.
  * General processing nodes which do not change the type of the moments
    may "inherit" the features from the parent node. For instance, slicing operator.
  * Conversions can be done easily in both of the above cases if the
    message conversion is defined in the moments class. For instance,
    GaussianMarkovChain to Gaussian and GaussianDriftingMarkovChain to
    Gaussian.

The base class will handle the necessary relevant tasks now so it is not
necessary to define the class for nodes. However, in the future some
calculations might be moved from nodes to statistics classes.

Currently, only GaussianStatistics and GaussianMarkovChainStatistics are
implemented.
ENH: Add shape check for exponential family node phi update method

This check will raise an error if there's something wrong with the
shape.
ENH: Add initialization by mean and covariance to Gaussian array node
DEMO: Improve drifting LSSM demos

Improve the initializations and the iteration to be a bit more robust.
Utilize the plot monitoring of the nodes so the user can see the
distributions during the iteration.
BUG: Fix sharex handling in timeseries plotting

Previously, if plotting multiple sets of timeseries, the axes were
cleared before each plot. Thus only one set of timeseries was shown, the
last one. This was because, if subplot gets a different value for sharex
than before, it will create a new axes object. Now, the plotting is done
so that the sharex value remains the same and the user can thus plot
multiple sets of timeseries in one figure.