IMP logo
IMP Reference Guide  2.6.0
The Integrative Modeling Platform
vonMisesSufficient.h
Go to the documentation of this file.
1 /**
2  * \file IMP/isd/vonMisesSufficient.h \brief Normal distribution of Function
3  *
4  * Copyright 2007-2016 IMP Inventors. All rights reserved.
5  */
6 
7 #ifndef IMPISD_VON_MISES_SUFFICIENT_H
8 #define IMPISD_VON_MISES_SUFFICIENT_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 //! vonMisesSufficient
20 /** Probability density function and -log(p) of von Mises distribution
21  of N iid von Mises observations, provided through their sufficient
22  statistics.
23  This is much more efficient than multiplying N von Mises densities.
24  \f[ f(\chi|N, R_0, \chi_{exp}, \kappa) =
25  \frac{\exp \left(R_0 \kappa \cos (\chi - \chi_{exp})\right)}
26  {2\pi I_0(\kappa)^N} \f]
27  where
28  \f[ R = \sqrt{\left(\sum_{i=1}^N \cos \chi_{exp}^i\right)^2
29  + \left(\sum_{i=1}^N \cos \chi_{exp}^i\right)^2} \f]
30  \f[ \exp (i \chi_{exp}) = \frac{1}{R} \sum_{j=1}^N \exp(i \chi_{exp}^j) \f]
31  If \f$N=1\f$ and \f$\mu_1=\mu_2\f$ this reduces to the original von Mises
32  distribution with known mean and concentration.
33  \note derivative with respect to the mean \f$\chi_{exp}\f$ is not provided.
34  */
35 
36 class vonMisesSufficient : public Object {
37  public:
38  /** compute von Mises given the sufficient statistics
39  \param[in] chi
40  \param[in] N number of observations
41  \param[in] R0 component of N observations on the x axis (\f$R_0\f$)
42  \param[in] chiexp mean (\f$\chi_{exp}\f$)
43  \param[in] kappa concentration
44  */
45  vonMisesSufficient(double chi, unsigned N, double R0, double chiexp,
46  double kappa)
47  : Object("von Mises sufficient %1%"),
48  x_(chi),
49  R0_(R0),
50  chiexp_(chiexp) {
51  N_ = N;
52  force_set_kappa(kappa);
53  }
54 
55  /** compute von Mises given the raw observations
56  * this is equivalent to calling get_sufficient_statistics and then the other
57  * constructor.
58  \param[in] chi
59  \param[in] obs a list of observed angles (in radians).
60  \param[in] kappa concentration
61  */
62  vonMisesSufficient(double chi, Floats obs, double kappa)
63  : Object("von Mises sufficient %1%"), x_(chi) {
64  Floats stats = get_sufficient_statistics(obs);
65  N_ = (unsigned)stats[0];
66  R0_ = stats[1];
67  chiexp_ = stats[2];
68  force_set_kappa(kappa);
69  }
70 
71  /* energy (score) functions, aka -log(p) */
72  virtual double evaluate() const {
73  return logterm_ - R0_ * kappa_ * cos(x_ - chiexp_);
74  }
75 
76  virtual double evaluate_derivative_x() const {
77  return R0_ * kappa_ * sin(x_ - chiexp_);
78  }
79 
80  virtual double evaluate_derivative_kappa() const {
81  return -R0_ * cos(x_ - chiexp_) + double(N_) * I1_ / I0_;
82  }
83 
84  /* probability density function */
85  virtual double density() const {
86  return exp(R0_ * kappa_ * cos(x_ - chiexp_)) / (2 * IMP::PI * I0N_);
87  }
88 
89  /* getting parameters */
90  double get_x() { return x_; }
91  double get_R0() { return R0_; }
92  double get_chiexp() { return chiexp_; }
93  double get_N() { return N_; }
94  double get_kappa() { return kappa_; }
95 
96  /* change of parameters */
97  void set_x(double x) { x_ = x; }
98 
99  void set_R0(double R0) { R0_ = R0; }
100 
101  void set_chiexp(double chiexp) { chiexp_ = chiexp; }
102 
103  void set_N(unsigned N) {
104  N_ = N;
105  I0N_ = pow(I0_, static_cast<int>(N_));
106  logterm_ = log(2 * IMP::PI * I0N_);
107  }
108 
109  void set_kappa(double kappa) {
110  if (kappa_ != kappa) {
111  force_set_kappa(kappa);
112  }
113  }
114 
115  //! Compute sufficient statistics from a list of observations.
116  /** See Mardia and El-Atoum, "Bayesian inference for the von Mises-Fisher
117  distribution ", Biometrika, 1967.
118  \return the number of observations, \f$R_0\f$ (the component on the
119  x axis) and \f$\chi_{exp}\f$
120  */
122  unsigned N = data.size();
123  // mean cosine
124  double cosbar = 0;
125  double sinbar = 0;
126  for (unsigned i = 0; i < N; ++i) {
127  cosbar += cos(data[i]);
128  sinbar += sin(data[i]);
129  }
130  double R = sqrt(cosbar * cosbar + sinbar * sinbar);
131  double chi = acos(cosbar / R);
132  if (sinbar < 0) chi = -chi;
133  Floats retval(3);
134  retval[0] = N;
135  retval[1] = R;
136  retval[2] = chi;
137  return retval;
138  }
139 
141  /*IMP_OBJECT_INLINE(vonMisesSufficient, out << "vonMisesSufficient: " << x_
142  << ", " << N_ << ", " << R0_ << ", " << chiexp_
143  << ", " << kappa_ <<std::endl, {});*/
144 
145  private:
146  double x_, R0_, chiexp_, kappa_, I0_, I1_, logterm_, I0N_;
147  unsigned N_;
148  void force_set_kappa(double kappa) {
149  kappa_ = kappa;
150  I0_ = double(boost::math::cyl_bessel_i(0, kappa));
151  I1_ = double(boost::math::cyl_bessel_i(1, kappa));
152  I0N_ = pow(I0_, static_cast<int>(N_));
153  logterm_ = log(2 * IMP::PI * I0N_);
154  }
155 };
156 
157 IMPISD_END_NAMESPACE
158 
159 #endif /* IMPISD_VON_MISES_SUFFICIENT_H */
Declare an efficient stl-compatible map.
static const double PI
the constant pi
Various useful constants.
#define IMP_OBJECT_METHODS(Name)
Define the basic things needed by any Object.
Definition: object_macros.h:25
static Floats get_sufficient_statistics(Floats data)
Compute sufficient statistics from a list of observations.
vonMisesSufficient(double chi, unsigned N, double R0, double chiexp, double kappa)
Storage of a model, its restraints, constraints and particles.
Various general useful macros for IMP.
Common base class for heavy weight IMP objects.
Definition: Object.h:106
vonMisesSufficient(double chi, Floats obs, double kappa)