00001 // This is mul/pdf1d/pdf1d_exponential_sampler.cxx 00002 00003 //: 00004 // \file 00005 // \author Tim Cootes 00006 // \brief Sampler class for Univariate exponential distributions 00007 00008 #include "pdf1d_exponential_sampler.h" 00009 00010 #include <vcl_cassert.h> 00011 #include <vcl_cmath.h> 00012 #include <pdf1d/pdf1d_exponential.h> 00013 00014 //======================================================================= 00015 // Dflt ctor 00016 //======================================================================= 00017 00018 pdf1d_exponential_sampler::pdf1d_exponential_sampler(): 00019 rng_(9667566ul) 00020 { 00021 } 00022 00023 //======================================================================= 00024 // Destructor 00025 //======================================================================= 00026 00027 pdf1d_exponential_sampler::~pdf1d_exponential_sampler() 00028 { 00029 } 00030 00031 00032 //======================================================================= 00033 // Method: is_a 00034 //======================================================================= 00035 00036 vcl_string pdf1d_exponential_sampler::is_a() const 00037 { 00038 return vcl_string("pdf1d_exponential_sampler"); 00039 } 00040 00041 //======================================================================= 00042 // Method: is_class 00043 //======================================================================= 00044 00045 bool pdf1d_exponential_sampler::is_class(vcl_string const& s) const 00046 { 00047 return pdf1d_sampler::is_class(s) || s==pdf1d_exponential_sampler::is_a(); 00048 } 00049 00050 //======================================================================= 00051 // Method: clone 00052 //======================================================================= 00053 00054 pdf1d_sampler* pdf1d_exponential_sampler::clone() const 00055 { 00056 return new pdf1d_exponential_sampler(*this); 00057 } 00058 //======================================================================= 00059 00060 void pdf1d_exponential_sampler::reseed(unsigned long seed) 00061 { 00062 rng_.reseed(seed); 00063 } 00064 //======================================================================= 00065 00066 //: Set model for which this is an instance 00067 // Error check that it is an axis exponential. 00068 void pdf1d_exponential_sampler::set_model(const pdf1d_pdf& model) 00069 { 00070 assert(model.is_class("pdf1d_exponential")); 00071 // cannot use dynamic_cast<> without rtti - PVr 00072 // rtti currently turned off 00073 pdf1d_sampler::set_model(model); 00074 } 00075 00076 //======================================================================= 00077 00078 double pdf1d_exponential_sampler::sample() 00079 { 00080 const pdf1d_exponential & exponential = static_cast<const pdf1d_exponential &>(model()); 00081 double L = exponential.lambda(); 00082 00083 return -1.0*vcl_log(rng_.drand64(0,1))/L; 00084 } 00085 00086 //: Fill x with samples possibly chosen so as to represent the distribution 00087 // 5 or fewer samples requested, they are spaced out equally. 00088 void pdf1d_exponential_sampler::regular_samples(vnl_vector<double>& x) 00089 { 00090 int n = x.size(); 00091 const pdf1d_exponential & exponential = static_cast<const pdf1d_exponential &>(model()); 00092 double L = exponential.lambda(); 00093 00094 // CDF = 1-exp(-Lx) 00095 // Require CDF(x) = (1+i)/(n+1) 00096 // Thus exp(-Lx) = (n-i)/(n+1) 00097 for (int i=0;i<n;++i) 00098 x[i] = -1.0*vcl_log(double(n-i)/(n+1))/L; 00099 } 00100 00101 00102 //======================================================================= 00103 00104 00105 //: Return a reference to the pdf model 00106 // This is properly cast. 00107 const pdf1d_exponential& pdf1d_exponential_sampler::exponential() const 00108 { 00109 return static_cast<const pdf1d_exponential&>(model()); 00110 } 00111