stochasticadaptiverungekutta.hpp

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

#ifndef INDII_ML_SDE_STOCHASTICADAPTIVERUNGEKUTTA_HPP
#define INDII_ML_SDE_STOCHASTICADAPTIVERUNGEKUTTA_HPP

#include "StochasticDifferentialModel.hpp"
#include "../ode/NumericalSolver.hpp"
#include "../aux/vector.hpp"
#include "../aux/matrix.hpp"
#include "../aux/WienerProcess.hpp"

#include <gsl/gsl_odeiv.h>
#include <map>

namespace indii {
  namespace ml {
    namespace sde {
/**
 * Stochastic Adaptive Runge-Kutta method for solving a system of
 * stochastic differential equations.
 *
 * @author Lawrence Murray <[e-mail]>
 * @version $Rev: 404 $
 * @date $Date: 2008-03-05 14:52:55 +0000 (Wed, 05 Mar 2008) $
 *
 * This class numerically solves StochasticDifferentialModel models
 * defining a system of stochastic differential equations using an
 * adaptive time step 8th/9th order (4th/4.5th order for stochastic
 * systems) Runge-Kutta Prince-Dormand method, as implemented in the
 * @ref GSL "GSL".
 *
 * The general usage idiom is as for
 * indii::ml::ode::AdaptiveRungeKutta. The class will automatically
 * keep track of Wiener process increments.
 *
 * @section StochasticAdaptiveRungeKutta_references References
 *
 * @anchor Wilkie2004
 * Wilkie, J. Numerical methods for stochastic differential
 * equations. <i>Physical Review E</i>, <b>2004</b>, <i>70</i>
 *
 * @anchor Sarkka2006
 * Särkkä, S. Recursive Bayesian Inference on Stochastic Differential
 * Equations. PhD thesis, <i>Helsinki University of Technology</i>,
 * <b>2006</b>.
 */
class StochasticAdaptiveRungeKutta :
    public indii::ml::ode::NumericalSolver {
public:
  /**
   * Constructor.
   *
   * @param model Model to estimate.
   *
   * The time is initialised to zero, but the state is uninitialised
   * and should be set with setVariable() or setState().
   */
  StochasticAdaptiveRungeKutta(StochasticDifferentialModel* model);

  /**
   * Constructor.
   *
   * @param model Model to estimate.
   * @param y0 Initial state.
   *
   * The time is initialised to zero and the state to that given.
   */
  StochasticAdaptiveRungeKutta(StochasticDifferentialModel* model,
      const indii::ml::aux::vector& y0);

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

  virtual double step(double upper);

  virtual double stepBack(double lower);

protected:
  virtual void reset();

private:
  /**
   * Model.
   */
  StochasticDifferentialModel* model;

  /**
   * Wiener process system noise.
   */
  indii::ml::aux::WienerProcess<double> W;

  /**
   * Wiener process increments at future times,
   * \f$\Delta\mathbf{W}(t_1),\Delta\mathbf{W}(t_2),\ldots\f$ indexed
   * by \f$t_1,t_2,\ldots > t\f$.
   */
  std::map<double,indii::ml::aux::vector> noisePath;

  /**
   * Last Wiener process increment sampled by sampleNoise(). This is
   * stored for reasons of efficiency, such that it needn't be
   * retrieved in O(log n) from noisePath.
   */
  indii::ml::aux::vector deltaW;

  /**
   * GSL ordinary differential equations step type.
   */
  static const gsl_odeiv_step_type* gslStepType;

  /**
   * Function for calculating derivatives when moving forward in time.
   */
  static indii::ml::ode::df_t gslForwardFunction;

  /**
   * Function for calculating derivatives when moving backward in
   * time.
   */
  static indii::ml::ode::df_t gslBackwardFunction;

  /**
   * Sample Wiener process noise. This conditions on previous samples
   * of the Wiener trajectory at future times resulting from time
   * steps rejected by the adaptive step size control. This ensures
   * that steps are rejected only due to error bounds being exceeded
   * and not due to an improbable but perfectly valid Wiener
   * trajectory.
   *
   * @param ts \f$t_s\f$; the time at which to sample the noise.
   *
   * Let \f$t\f$ be the current time and \f$t_s\f$ the time at which
   * to sample the noise.
   *
   * Of the stored Wiener increments at times \f$t_1,t_2,\ldots >
   * t\f$, let \f$t_l\f$ be the greatest time less than or equal to
   * \f$t_s\f$, and \f$t_u\f$ the least time greater than
   * \f$t_s\f$. \f$\Delta\mathbf{W}(t_l)\f$ and
   * \f$\Delta\mathbf{W}(t_u)\f$ are the corresponding Wiener
   * increments.
   *
   * If \f$t_l\f$ does not exist, set \f$t_l \leftarrow t\f$.
   *
   * Then, if \f$t_s \neq t_l\f$, draw \f$\mathbf{\delta} \sim
   * \mathcal{N}(0,\sqrt{t_s-t_l})\f$. If \f$t_u\f$ does not exist,
   * set \f$\Delta\mathbf{W}(t_s) \leftarrow
   * \mathbf{\delta}\f$. Otherwise, set:
   *
   * \f{eqnarray*}
   * \Delta\mathbf{W}(t_s) &\leftarrow& \frac{t_s - t_l}{t_u -
   * t_l}\Delta\mathbf{W}(t_u) + \mathbf{\delta} \\
   * \Delta\mathbf{W}(t_u) &\leftarrow& \frac{t_u - t_s}{t_u -
   * t_l}\Delta\mathbf{W}(t_u) - \mathbf{\delta} \\
   * \f}
   *
   * and store.
   *
   * @return \f$\Delta\mathbf{W}(t_s)\f$; increment of the noise at
   * the given time.
   */
  indii::ml::aux::vector sampleNoise(const double ts);

};

    }
  }
}

#endif