KJB
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
prob_pdf.h
Go to the documentation of this file.
1 /* ======================================================================== *
2  | |
3  | Copyright (c) 2007-2010, by members of University of Arizona Computer |
4  | Vision group (the authors) including |
5  | Kobus Barnard, Kyle Simek, Ernesto Brau. |
6  | |
7  | For use outside the University of Arizona Computer Vision group please |
8  | contact Kobus Barnard. |
9  | |
10  * ======================================================================== */
11 
12 /* $Id: prob_pdf.h 17970 2014-10-28 15:07:45Z ernesto $ */
13 
14 #ifndef PROB_PDF_H_INCLUDED
15 #define PROB_PDF_H_INCLUDED
16 
31 #include "m_cpp/m_special.h"
33 
34 namespace kjb {
35 
36 /*============================================================================
37  pdf functions; notice that we import all boost pdfs here as well
38  ----------------------------------------------------------------------------*/
39 
40 // importing the boost::math::pdf into this namespace so that we
41 // can call the boost ones using the kjb:: prefix.
42 using boost::math::pdf;
43 
44 #if BOOST_VERSION < 104600
45 
48 inline double pdf(const Geometric_distribution& dist, double x)
49 {
50  return boost::math::pdf<boost::math::negative_binomial, double>(dist, x);
51 }
52 #endif
53 
54 
58 inline double pdf(const Log_normal_distribution& dist, double x)
59 {
60  return pdf(dist.dist_, log(x)) / x;
61 }
62 
66 template<class T>
67 inline
68 double pdf(const Categorical_distribution<T>& dist, const T& x)
69 {
70  typename std::map<T, double>::const_iterator i = dist.db.find(x);
71  return i == dist.db.end() ? 0.0 : i->second/dist.total_weight_;
72 }
73 
77 double pdf(const MV_gaussian_distribution& P, const Vector& x);
78 
82 template<class Distribution>
83 double pdf
84 (
86  const typename
88 )
89 {
90  double f = 0.0;
91  for(size_t i = 0; i < dist.dists.size(); i++)
92  {
93  f += (pdf(dist.coeffs, i + 1) * pdf(dist.dists[i], x));
94  }
95 
96  return f;
97 }
98 
99 
103 inline double pdf(const Gaussian_distribution& P, double x)
104 {
105  return boost::math::pdf(P, x);
106 }
107 
108 // forward declaration
109 inline double log_pdf(const Beta_binomial_distribution& dist, double k);
110 
112 inline double pdf(const Beta_binomial_distribution& dist, double k)
113 {
114  return exp(log_pdf(dist, k));
115 }
116 
117 // forward declaration
118 template <size_t D>
119 inline double log_pdf(const Von_mises_fisher_distribution<D>& dist, const kjb::Vector_d<D>& v);
120 
122 template <size_t D>
123 inline
125 {
126  return exp(log_pdf(dist, v));
127 }
128 
129 // forward declaration
130 template <size_t D>
131 inline double log_pdf(const Uniform_sphere_distribution<D>& P, const kjb::Vector_d<D>& x);
132 
134 template <size_t D>
135 inline double pdf(const Uniform_sphere_distribution<D>& P, const kjb::Vector_d<D>& x)
136 {
137  return exp(log_pdf(P,x));
138 }
139 
140 /*============================================================================
141  cdf functions; notice that we import all boost cdfs here as well
142  ----------------------------------------------------------------------------*/
143 
144 // importing the boost::math::cdf into this namespace so that we
145 // can call the boost ones using the kjb:: prefix.
146 using boost::math::cdf;
147 
151 inline double cdf(const Log_normal_distribution& dist, double x)
152 {
153  return cdf(dist.dist_, log(x));
154 }
155 
159 template<class T>
160 double cdf(const Categorical_distribution<T>& dist, const T& x)
161 {
162  // TODO: compute this in O(log n) time using the precomputed cdf_ vector and std::upper_bound. This code should just work, but no time to test right now:
163  // Update: actually, this probably won't work, because the CDF isn't guaranteed to be stored in key order. this needs some thought, but getting a proper O(log n) algorithm shouldn't be that hard... -ksimek, Aug 13, 2013
164 // typedef Categorical_distribution<T> Cd;
165 // return *--std::upper_bound(
166 // dist.cdf_.begin(),
167 // dist.cdf_.end(),
168 // x,
169 // Cd::cdf_key_less_than_scalar());
170 
171  // compute CDF on-the-fly in O(n) time.
172  double F = 0.0;
173  typename std::map<T, double>::const_iterator y = dist.db.upper_bound(x);
174  for(typename std::map<T, double>::const_iterator i = dist.db.begin();
175  i != y; i++)
176  {
177  F += i->second / dist.total_weight_;
178  }
179 
180  return F;
181 }
182 
186 template<class Distribution>
187 double cdf
188 (
190  const typename
192 )
193 {
194  double f = 0.0;
195  for(size_t i = 0; i < dist.dists.size(); i++)
196  {
197  f += (cdf(dist.coeffs, i + 1) * cdf(dist.dists[i], x));
198  }
199 
200  return f;
201 }
202 
203 /*============================================================================
204  quantile functions; notice that we import all boost quantile here as well
205  ----------------------------------------------------------------------------*/
206 
207 using boost::math::quantile;
208 
209 /*============================================================================
210  log-pdf functions; notice the generic one...
211  ----------------------------------------------------------------------------*/
212 
219 template<class Distribution>
220 inline double log_pdf
221 (
222  const Distribution& P,
224 )
225 {
226  return std::log(pdf(P, x));
227 }
228 
229 
233 inline double log_pdf(const Beta_binomial_distribution& dist, double k)
234 {
235  const double lcoeff = log_binomial_coefficient(dist.n, k);
236  const double na = k + dist.alpha;
237  const double nb = dist.n - k + dist.beta;
238 
239  const double lnumer = lgamma(na) + lgamma(nb) - lgamma(na+nb);
240  const double ldenom = lgamma(dist.alpha) + lgamma(dist.beta) - lgamma(dist.alpha + dist.beta);
241 
242  return lcoeff + lnumer - ldenom;
243 }
244 
251 inline double log_pdf(const Gaussian_distribution& P, double x)
252 {
253  double mu = P.mean();
254  double variance = P.standard_deviation() * P.standard_deviation();
255 
256  return -(0.5 * log(2 * M_PI * variance))
257  - (((x - mu) * (x - mu)) / (2 * variance));
258 }
259 
266 inline double log_pdf(const Laplace_distribution& P, double x)
267 {
268  double mu = P.location();
269  double b = P.scale();
270 
271  return -log(2 * b) - (fabs(x - mu) / b);
272 }
273 
280 inline double log_pdf(const Binomial_distribution& P, double k)
281 {
282  const double p = P.success_fraction();
283  const double n = P.trials();
284  double result = log_binomial_coefficient(n, k);
285  result += k * log(p) + (n - k) * log(1-p);
286  return result;
287 }
288 
295 double log_pdf(const MV_gaussian_distribution& P, const Vector& x);
296 
297 
304 template <size_t D>
305 inline double log_pdf(const Uniform_sphere_distribution<D>& /*dist*/, const kjb::Vector_d<D>& /*v*/)
306 {
307  return -log(kjb_c::unit_sphere_surface_area(D));
308 }
309 
310 
311 template <size_t D>
313 {
314  static const double LOG_2_PI = log(2*M_PI);
315 
316  double kappa_ = dist.kappa();
317  const kjb::Vector_d<D>& mu_ = dist.mu();
318 
319  double log_exponential = kappa_ * dot(v.normalized(), mu_);
320  if(D == 3)
321  {
322  //TODO: compute normalization inside constructor
323  double log_normalization = log(kappa_) - (LOG_2_PI + log(exp(kappa_) - exp(-kappa_)));
324 
325 #ifdef TEST
326  // math check - implement normalization in 3 other ways
327  double test_1 = log(kappa_ / (2 * M_PI * (exp(kappa_) - exp(-kappa_))));
328  double test_2 = log(kappa_ / (4* M_PI * std::sinh(kappa_)));
329  // general form
330  double test_3 = (D/2.0 - 1) * log(kappa_) - (D/2.0)*log(2 * M_PI) - log(boost::math::cyl_bessel_i(D/2.0-1, kappa_));
331 
332  // THESE ALL MATCH AS OF May 18, 2012
333  assert(fabs(log_normalization - test_1) < FLT_EPSILON);
334  assert(fabs(log_normalization - test_2) < FLT_EPSILON);
335  assert(fabs(log_normalization - test_3) < FLT_EPSILON);
336 #endif
337  return log_exponential + log_normalization;
338  }
339  else
340  {
341  //KJB_THROW(kjb::Not_implemented);
342 
343  // THIS SHOULD BE RIGHT AND IS TESTED FOR THE D==3 case:
344  double log_normalization = (D/2.0 - 1) * log(kappa_) - (D/2.0)*log(2 * M_PI) - log(boost::math::cyl_bessel_i(D/2.0-1, kappa_));
345  return log_exponential + log_normalization;
346  // implementing this will be easy once the 3D version is confirmed:
347  // just use the general form of the normalization factor (test_3 above)
348  // and use the same equation for the log_exponential
349  }
350 }
351 
353 double pdf
354 (
355  const Chinese_restaurant_process& cpr,
357 );
358 
360 double log_pdf
361 (
362  const Chinese_restaurant_process& cpr,
364 );
365 
366 } //namespace kjb
367 
368 #endif /*PROB_PDF_H_INCLUDED */
369 
const kjb::Vector_d< D > & mu() const
Definition: prob_distribution.h:996
Definition of various standard probability distributions.
double log_binomial_coefficient(size_t n, size_t k)
Definition: m_special.cpp:25
double kappa() const
Definition: prob_distribution.h:998
std::vector< std::vector< size_t > > Type
Definition: prob_distribution.h:1019
This class implements a mixture distribution. In other words, it is the sum of a finite number of fra...
Definition: prob_distribution.h:867
for k
Definition: APPgetLargeConnectedEdges.m:61
Log-normal distribution.
Definition: prob_distribution.h:827
double cdf(const Mixture_distribution< Distribution > &dist, const typename Distribution_traits< Mixture_distribution< Distribution > >::type &x)
Computes the CDF of a mixture distribution at x.
Definition: prob_pdf.h:188
Vector_d< D > normalized() const
Non-mutating (functionally pure) version of normalize()
Definition: m_vector_d.impl.h:384
boost::math::normal Gaussian_distribution
Definition: prob_distribution.h:66
Definition: prob_distribution.h:984
double pdf(const MV_gaussian_distribution &P, const Vector &x)
Computes the joint PDF of a multivariate Gaussian at x.
Definition: prob_pdf.cpp:32
A categorical distribution.
Definition: prob_distribution.h:131
double beta
Definition: prob_distribution.h:96
boost::math::laplace Laplace_distribution
Definition: prob_distribution.h:67
double dist(const pt a, const pt b)
compute approx. Great Circle distance between two UTM points
Definition: layer.cpp:45
Definition: prob_distribution.h:966
double n
Definition: prob_distribution.h:96
double type
Definition: prob_distribution.h:111
x
Definition: APPgetLargeConnectedEdges.m:100
Generic traits struct for all distributions.
Definition: prob_distribution.h:109
double log_pdf(const MV_gaussian_distribution &P, const Vector &x)
Computes the log PDF a multivariate normal distribution at x.
Definition: prob_pdf.cpp:64
double pdf(const Uniform_sphere_distribution< D > &P, const kjb::Vector_d< D > &x)
Evaluate the probability density function of a uniform sphere distribution (x's magnitude is ignored ...
Definition: prob_pdf.h:135
Definition: prob_distribution.h:79
boost::math::binomial Binomial_distribution
Definition: prob_distribution.h:62
#define M_PI
Definition: fft.cpp:206
double cdf(const Log_normal_distribution &dist, double x)
Computes the cdf of a log-normal distribution at x.
Definition: prob_pdf.h:151
Definition: g_quaternion.h:37
get the indices of edges in each direction for i
Definition: APPgetLargeConnectedEdges.m:48
long int dot(const Int_vector &op1, const Int_vector &op2)
Returns dot product of this and op2.
Definition: l_int_vector.h:1532
D
Definition: APPgetLargeConnectedEdges.m:106
Definition: prob_distribution.h:89
double alpha
Definition: prob_distribution.h:96