IMP  2.0.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-2013 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 
37 {
38  public:
39  /** compute von Mises given the sufficient statistics
40  \param[in] chi
41  \param[in] N number of observations
42  \param[in] R0 component of N observations on the x axis (\f$R_0\f$)
43  \param[in] chiexp mean (\f$\chi_{exp}\f$)
44  \param[in] kappa concentration
45  */
46  vonMisesSufficient(double chi, unsigned N, double R0, double chiexp,
47  double kappa):
48  base::Object("von Mises sufficient %1%"), x_(chi), R0_(R0), chiexp_(chiexp)
49 
50  {
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  base::Object("von Mises sufficient %1%"), x_(chi)
64  {
65  Floats stats = get_sufficient_statistics(obs);
66  N_= (unsigned) stats[0];
67  R0_ = stats[1];
68  chiexp_ = stats[2];
69  force_set_kappa(kappa);
70  }
71 
72  /* energy (score) functions, aka -log(p) */
73  virtual double evaluate() const
74  {
75  return logterm_ - R0_*kappa_*cos(x_-chiexp_);
76  }
77 
78  virtual double evaluate_derivative_x() const
79  {
80  return R0_*kappa_*sin(x_-chiexp_) ;
81  }
82 
83  virtual double evaluate_derivative_kappa() const
84  {
85  return - R0_ * cos(x_-chiexp_) + double(N_) * I1_/I0_ ;
86  }
87 
88  /* probability density function */
89  virtual double density() const
90  {
91  return exp(R0_*kappa_*cos(x_-chiexp_))/(2*IMP::PI*I0N_);
92  }
93 
94  /* getting parameters */
95  double get_x() { return x_; }
96  double get_R0() { return R0_; }
97  double get_chiexp() { return chiexp_; }
98  double get_N() { return N_; }
99  double get_kappa() { return kappa_; }
100 
101  /* change of parameters */
102  void set_x(double x) {
103  x_=x;
104  }
105 
106  void set_R0(double R0) {
107  R0_=R0;
108  }
109 
110  void set_chiexp(double chiexp) {
111  chiexp_=chiexp;
112  }
113 
114  void set_N(unsigned N){
115  N_=N;
116  I0N_=pow(I0_, static_cast<int>(N_));
117  logterm_ = log(2*IMP::PI*I0N_);
118  }
119 
120  void set_kappa(double kappa) {
121  if (kappa_ != kappa) {
122  force_set_kappa(kappa);
123  }
124  }
125 
126  //! Compute sufficient statistics from a list of observations.
127  /** See Mardia and El-Atoum, "Bayesian inference for the von Mises-Fisher
128  distribution ", Biometrika, 1967.
129  \return the number of observations, \f$R_0\f$ (the component on the
130  x axis) and \f$\chi_{exp}\f$
131  */
132  static Floats get_sufficient_statistics(Floats data)
133  {
134  unsigned N = data.size();
135  //mean cosine
136  double cosbar=0;
137  double sinbar=0;
138  for (unsigned i=0; i<N; ++i){
139  cosbar += cos(data[i]);
140  sinbar += sin(data[i]);
141  }
142  double R=sqrt(cosbar*cosbar + sinbar*sinbar);
143  double chi=acos(cosbar/R);
144  if (sinbar < 0) chi=-chi;
145  Floats retval (3);
146  retval[0]=N;
147  retval[1]=R;
148  retval[2]=chi;
149  return retval;
150  }
151 
152  IMP_OBJECT_INLINE(vonMisesSufficient, out << "vonMisesSufficient: " << x_
153  << ", " << N_ << ", " << R0_ << ", " << chiexp_
154  << ", " << kappa_ <<std::endl, {});
155 
156  private:
157  double x_,R0_,chiexp_,kappa_,I0_,I1_,logterm_,I0N_;
158  unsigned N_;
159  void force_set_kappa(double kappa) {
160  kappa_ = kappa;
161  I0_ = double(boost::math::cyl_bessel_i(0, kappa));
162  I1_ = double(boost::math::cyl_bessel_i(1, kappa));
163  I0N_=pow(I0_, static_cast<int>(N_));
164  logterm_ = log(2*IMP::PI*I0N_);
165  }
166 };
167 
168 IMPISD_END_NAMESPACE
169 
170 #endif /* IMPISD_VON_MISES_SUFFICIENT_H */