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_