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