core/vpdl/vpdt/vpdt_log_probability.h
Go to the documentation of this file.
00001 // This is core/vpdl/vpdt/vpdt_log_probability.h
00002 #ifndef vpdt_log_probability_h_
00003 #define vpdt_log_probability_h_
00004 //:
00005 // \file
00006 // \author Matthew Leotta
00007 // \brief The basic functions for log of probability calculation
00008 // \date March 13, 2009
00009 //
00010 // These functions provide default implementations for various log of
00011 // probability calculation functions.  They are written in terms of distribution
00012 // member functions
00013 //
00014 // \verbatim
00015 //  Modifications
00016 //   None
00017 // \endverbatim
00018 
00019 #include <vpdl/vpdt/vpdt_dist_traits.h>
00020 #include <vnl/vnl_math.h>
00021 #include <vcl_limits.h>
00022 
00023 //: Compute the log of the unnormalized density
00024 template <class dist>
00025 inline typename vpdt_dist_traits<dist>::scalar_type
00026 vpdt_log_density(const dist& d,
00027                  const typename vpdt_dist_traits<dist>::field_type& pt)
00028 {
00029   typedef typename vpdt_dist_traits<dist>::scalar_type T;
00030   T density = d.density(pt);
00031   if (density <= T(0))
00032     return vcl_numeric_limits<T>::infinity();
00033 
00034   return static_cast<T>(vcl_log(density));
00035 }
00036 
00037 
00038 //: Compute the log of the normalized probability density
00039 template <class dist>
00040 inline typename vpdt_dist_traits<dist>::scalar_type
00041 vpdt_log_prob_density(const dist& d,
00042                       const typename vpdt_dist_traits<dist>::field_type& pt)
00043 {
00044   typedef typename vpdt_dist_traits<dist>::scalar_type T;
00045   T norm = d.norm_const();
00046   if (vnl_math_isinf(norm))
00047     return -vcl_numeric_limits<T>::infinity();
00048 
00049   return static_cast<T>(vcl_log(norm) + vpdt_log_density(d,pt));
00050 }
00051 
00052 
00053 //: Compute the gradient of the log of the unnormalized density
00054 template <class dist>
00055 inline typename vpdt_dist_traits<dist>::scalar_type
00056 vpdt_gradient_log_density(const dist& d,
00057                           const typename vpdt_dist_traits<dist>::field_type& pt,
00058                           const typename vpdt_dist_traits<dist>::vector_type& g)
00059 {
00060   typedef typename vpdt_dist_traits<dist>::scalar_type T;
00061   T density = d.gradient_density(pt,g);
00062   if (density <= T(0)) {
00063     vpdt_fill(g,T(0));
00064     return vcl_numeric_limits<T>::infinity();
00065   }
00066 
00067   g /= density;
00068   return static_cast<T>(vcl_log(density));
00069 }
00070 
00071 
00072 #endif // vpdt_log_probability_h_