wienerprocess.hpp

dysii是一款非常出色的滤波函数库

#ifndef INDII_ML_AUX_WIENERPROCESS_HPP
#define INDII_ML_AUX_WIENERPROCESS_HPP

#include "StochasticProcess.hpp"

namespace indii {
  namespace ml {
    namespace aux {

/**
 * Multivariate Wiener process.
 *
 * @author Lawrence Murray <[e-mail]>
 * @version $Rev: 366 $
 * @date $Date: 2008-01-22 17:40:25 +0000 (Tue, 22 Jan 2008) $
 *
 * The Wiener process, \f$W(t)\f$, is a diffusion process where
 * \f$W(t+\Delta t) - W(t) \sim \mathcal{N}(0, \sqrt{\Delta t})\f$, with
 * the properties of:
 *
 * @li being Markovian,
 * @li having continuous sample paths which are nowhere
 * differentiable, and
 * @li having highly irregular sample paths, with mean approaching
 * zero but variance approaching infinity as \f$\Delta t \rightarrow
 * \infty\f$.
 * 
 * Given the properties of deterministic functions, point 2 seems
 * self-contradicting, but such is the nature of stochastic processes!
 *
 * This implementation is generalised to allow a nonzero drift
 * \f$\mathbf{\mu}\f$ (expectation per unit time), and arbitrary
 * diffusion \f$\Sigma\f$ (covariance per unit time), so that the
 * process is given by:
 *
 * \f[
 * d\mathbf{x} = \mathbf{\mu}\,dt + L\,d\mathbf{W},
 * \f]
 *
 * where \f$\Sigma = LL^T\f$ (Cholesky decomposition) and the
 * \f$dW_i\f$ are independent, standard (zero drift, unit variance per
 * unit time) Wiener processes.
 */
template <class T = unsigned int>
class WienerProcess : public StochasticProcess<T> {
public:
  /**
   * Default constructor.
   *
   * Initialises the Wiener process with zero dimensions. This should
   * generally only be used when the object is to be restored from a
   * serialization.
   */
  WienerProcess();

  /**
   * Construct Wiener process with given drift and diffusion.
   *
   * @param mu \f$\mathbf{\mu}\f$; drift of the process.
   * @param sigma \f$\Sigma\f$; diffusion of the process.
   */
  WienerProcess(const vector& mu, const symmetric_matrix& sigma);

  /**
   * Construct Wiener process of given dimensionality.
   *
   * @param N \f$N\f$; number of dimensions in the process.
   *
   * The drift is initialised to zero and the diffusion to the
   * identity matrix, giving the standard Wiener process. These may be
   * modified with setDrift() and setDiffusion().
   */
  WienerProcess(unsigned int N);

  /**
   * Destructor.
   */
  virtual ~WienerProcess();

  /**
   * Assignment operator.
   */
  WienerProcess<T>& operator=(const WienerProcess<T>& o);
  
  virtual unsigned int getDimensions() const;

  virtual void setDimensions(const unsigned int N,
      const bool preserve = false);

  /**
   * Get the drift of the process.
   *
   * @return The drift of the process.
   */
  virtual const vector& getDrift() const;

  /**
   * Get the diffusion of the process.
   *
   * @return The diffusion of the process.
   */
  virtual const symmetric_matrix& getDiffusion() const;

  virtual const vector& getDrift();

  virtual const symmetric_matrix& getDiffusion();

  /**
   * Set the drift of the process. The new drift vector must have the
   * same size as the previous one.
   *
   * @param mu \f$\mathbf{\mu}\f$; drift of the process.
   */
  virtual void setDrift(const vector& mu);

  /**
   * Set the diffusion of the Gaussian. The new diffusion matrix
   * must have the same size as the previous one.
   *
   * @param sigma \f$\Sigma\f$; diffusion of the process.
   */
  virtual void setDiffusion(const symmetric_matrix& sigma);

  virtual vector getExpectation(const T delta);

  virtual symmetric_matrix getCovariance(const T delta);

  /**
   * Get the expected value of the process after a given time has
   * elapsed.
   *
   * @param delta \f$\Delta t\f$; elapsed time.
   *
   * @return \f$\mathbf{\mu}(\Delta t)\f$; expected value of the process
   * after time \f$\Delta t\f$.
   */
  virtual vector getExpectation(const T delta) const;

  /**
   * Get the covariance of the process after a given time has elapsed.
   *
   * @param delta \f$\Delta t\f$; elapsed time.
   *
   * @return \f$\Sigma(\Delta t)\f$; covariance of the process after
   * time \f$\Delta t\f$.
   */
  virtual symmetric_matrix getCovariance(const T delta) const;

  /**
   * Sample from the process after a given time has elapsed.
   *
   * @param delta \f$\Delta t\f$; elapsed time.
   *
   * This is performed using the following process:
   *
   * @li Let \f$\mathbf{z}\f$ be a vector of \f$N\f$ independent
   * normal variates with mean zero and variance \f$\Delta t\f$.
   *
   * @li Then the sample is \f$\mathbf{x} = \Delta t \mathbf{\mu} +
   * L\mathbf{z}\f$, where \f$L\f$ is the Cholesky decomposition of
   * the diffusion.
   *
   * @return The sample \f$\mathbf{x}\f$
   */
  virtual vector sample(const T delta);

  /**
   * Calculate the density of the distribution at a given point after
   * a given time has elapsed:
   *
   * \f[p(\mathbf{x},\Delta t) =
   *   \frac{1}{\sqrt{(2\pi\Delta t)^N|\Sigma|}}
   *   \exp\big({-\frac{1}{2\Delta t}(\mathbf{x}-\mathbf{\mu})^T \Sigma^{-1} 
   *   (\mathbf{x}-\mathbf{\mu})\big)}
   * \f]
   *
   * @param delta \f$\Delta t\f$; elapsed time.
   * @param x \f$\mathbf{x}\f$; the point at which to calculate the
   * density.
   *
   * @return The density of the distribution at \f$\mathbf{x}\f$ after
   * time \f$\Delta t\f$.
   */
  virtual double densityAt(const T delta, const vector& x);

  /**
   * @deprecated Use densityAt() instead.
   */
  double calculateDensity(const T delta, const vector& x);

protected:
  /**
   * Called when changes are made to the process, such as when the
   * drift or diffusion is changed. This allows pre-calculations to be
   * refreshed.
   */
  virtual void dirty();

private:
  /**
   * \f$N\f$; number of dimensions.
   */
  unsigned int N;

  /**
   * \f$\mu\f$; drift.
   */
  vector mu;

  /**
   * \f$\Sigma\f$; diffusion.
   */
  symmetric_matrix sigma;

  /**
   * \f$L\f$ where \f$\Sigma = LL^T\f$ is the Cholesky decomposition
   * of the covariance matrix.
   */
  lower_triangular_matrix L;

  /**
   * \f$\Sigma^{-1}\f$
   */
  symmetric_matrix sigmaI;

  /**
   * Diagonal of \f$\Sigma^{-1}\f$.
   */
  vector sigmaIDiag;

  /**
   * \f$\det(\Sigma)\f$
   */
  double sigmaDet;

  /**
   * \f$\frac{1}{Z}\f$
   */
  double ZI;

  /**
   * Is the expectation zero?
   */
  bool isMuZero;

  /**
   * Is the diffusion matrix diagonal?
   */
  bool isSigmaIDiagonal;

  /**
   * Has the Cholesky decomposition of the covariance matrix been computed?
   */
  bool haveCholesky;

  /**
   * Has the inverse of the covariance matrix been computed?
   */
  bool haveInverse;

  /**
   * Has the determinant of the covariance matrix been computed?
   */
  bool haveDeterminant;

  /**
   * Has the normalising constant for the density function been
   * computed (up to \f$1/(\Delta t)^{N/2}\f$)?
   */
  bool haveZ;

  /**
   * Have optimisations for density calculations been determined?
   */
  bool haveDensityOptimisations;

  /**
   * Calculate the Cholesky decomposition of the covariance matrix.
   */
  void calculateCholesky();

  /**
   * Calculate the inverse of the diffusion matrix. This exploits the
   * calculation of the Cholesky decomposition
   */
  void calculateInverse();

  /**
   * Calculate the determinant of the diffusion matrix,
   * \f$|\Sigma|\f$. This exploits the calculation of the Cholesky
   * decomposition \f$L\f$. For any matrices \f$A\f$ and \f$B\f$,
   * \f$|AB| = |A||B|\f$. Given that \f$\Sigma = LL^T\f$, \f$|\Sigma|
   * = |L|^2\f$. The determinant of a triangular matrix is simply the
   * product of the elements on its diagonal, so \f$|L|\f$ is easy to
   * calculate.
   */
  void calculateDeterminant();

  /**
   * Calculate part of the normalising constant, \f$\frac{1}{Z}\f$ for
   * the density function.
   *
   * \f[\sqrt{(2\pi)^N|\Sigma|}\f]
   *
   * where the complete constant requires the additional
   * time-dependent factor:
   *
   * \f[Z = (\Delta t)^{N/2}\f]
   */
  void calculateZ();

  /**
   * Calculate optimisations for the calculation of densities.
   */
  void calculateDensityOptimisations();

  /**
   * Serialize.
   */
  template<class Archive>
  void save(Archive& ar, const unsigned int version) const;

  /**
   * Restore from serialization.
   */
  template<class Archive>
  void load(Archive& ar, const unsigned int version);

  /*
   * Boost.Serialization requirements.
   */
  BOOST_SERIALIZATION_SPLIT_MEMBER()
  friend class boost::serialization::access;

};

    }
  }
}

#include "Random.hpp"

#include "boost/numeric/bindings/traits/ublas_matrix.hpp"
#include "boost/numeric/bindings/traits/ublas_vector.hpp"
#include "boost/numeric/bindings/lapack/lapack.hpp"
#include "boost/serialization/base_object.hpp"

#include <assert.h>

using namespace indii::ml::aux;

namespace ublas = boost::numeric::ublas;
namespace lapack = boost::numeric::bindings::lapack;

template <class T>
WienerProcess<T>::WienerProcess() : N(0), mu(0), sigma(0), L(0,0), sigmaI(0),
    sigmaIDiag(0), haveCholesky(false), haveInverse(false),
    haveDeterminant(false), haveZ(false), haveDensityOptimisations(false) {
  //