IMP  2.0.1
The Integrative Modeling Platform
vonMises.h
Go to the documentation of this file.
1 /**
2  * \file IMP/isd/vonMises.h \brief Normal distribution of Function
3  *
4  * Copyright 2007-2013 IMP Inventors. All rights reserved.
5  */
6 
7 #ifndef IMPISD_VON_MISES_H
8 #define IMPISD_VON_MISES_H
9 
10 #include <IMP/isd/isd_config.h>
11 #include <IMP/macros.h>
12 #include <IMP/Model.h>
13 #include <IMP/constants.h>
14 #include <math.h>
15 #include <boost/math/special_functions/bessel.hpp>
16 
17 IMPISD_BEGIN_NAMESPACE
18 
19 //! vonMises
20 /** Probability density function and -log(p) of von Mises distribution
21  \f[ f(x|\mu,\kappa) = \frac{\exp \left(\kappa \cos (x-\mu)\right)}{2\pi
22  I_0(\kappa)} \f]
23  This is the proper treatment for a "normally distributed" angle.
24  When \f$\kappa\f$ becomes infinite, the distribution tends to a gaussian,
25  and \f$\kappa = 1/\sigma^2\f$.
26  */
27 
28 class vonMises : public base::Object
29 {
30 public:
31  vonMises(double x, double mu, double kappa)
32  : base::Object("von Mises %1%"), x_(x), mu_(mu) {
33  force_set_kappa(kappa);
34  }
35 
36  /* energy (score) functions, aka -log(p) */
37  virtual double evaluate() const
38  {
39  return logterm_ - kappa_*cos(x_-mu_);
40  }
41 
42  virtual double evaluate_derivative_x() const
43  { return kappa_*sin(x_-mu_); }
44 
45  virtual double evaluate_derivative_mu() const
46  { return -kappa_*sin(x_-mu_); }
47 
48  virtual double evaluate_derivative_kappa() const
49  { return -cos(x_-mu_) + I1_/I0_; }
50 
51  /* probability density function */
52  virtual double density() const
53  {
54  return exp(kappa_*cos(x_-mu_))/(2*IMP::PI*I0_);
55  }
56 
57  /* change of parameters */
58  void set_x(double x) {
59  x_=x;
60  }
61  void set_mu(double mu) {
62  mu_=mu;
63  }
64  void set_kappa(double kappa) {
65  if (kappa_ != kappa) {
66  force_set_kappa(kappa);
67  }
68  }
69 
70  IMP_OBJECT_INLINE(vonMises, out << "vonMises: " << x_ << ", " << mu_
71  << ", " << kappa_ <<std::endl, {});
72 
73  private:
74  void force_set_kappa(double kappa) {
75  kappa_ = kappa;
76  I0_ = boost::math::cyl_bessel_i(0, kappa);
77  I1_ = boost::math::cyl_bessel_i(1, kappa);
78  logterm_ = log(2*IMP::PI*I0_);
79  }
80  double x_,mu_,kappa_,I0_,I1_,logterm_;
81 };
82 
83 IMPISD_END_NAMESPACE
84 
85 #endif /* IMPISD_VON_MISES_H */