contrib/brl/bbas/bsta/algo/bsta_fit_weibull.cxx Source File

Go to the documentation of this file.
00001 // This is brl/bbas/bsta/algo/bsta_fit_weibull.cxx
00002 #ifndef bsta_fit_weibull_cxx_
00003 #define bsta_fit_weibull_cxx_
00004 //:
00005 // \file
00006 #include "bsta_fit_weibull.h"
00007 #include <vnl/vnl_gamma.h>
00008 #include <vcl_cassert.h>
00009 
00010 
00011 //this function computes the model value for (mean/standard_dev)^2 as
00012 //a function of Weibull parameter k.
00013 static double mean_sigma_sq(double k)
00014 {
00015   double m = vnl_gamma((1+k)/k);
00016   m *=m;
00017   double den = vnl_gamma((2+k)/k)-m;
00018   return m/den;
00019 }
00020 //The cost function is ((mean/sigma)^2(k) - sample_mean^2/sample_variance)^2
00021 double bsta_weibull_cost_function::f(vnl_vector<double> const& x)
00022 {
00023   double k = x[0];
00024   assert(k>0);
00025   double r = mean_/std_dev_;
00026   r *= r;
00027   double r_weib = mean_sigma_sq(k);
00028   double ret = r_weib-r;
00029   return ret*ret;
00030 }
00031 
00032 void bsta_weibull_cost_function::gradf(vnl_vector<double> const& x,
00033                                        vnl_vector<double>& gradient)
00034 {
00035   double k = x[0];
00036   assert(k>0);
00037   double m = vnl_gamma((1+k)/k);
00038   m*=m;
00039   double v = vnl_gamma((2+k)/k);
00040   double tn = vnl_digamma((2+k)/k)-vnl_digamma((1+k)/k);
00041   double neu = 2.0*m*v*tn;
00042   double den = k*k*(m-v)*(m-v);
00043   double r = mean_/std_dev_;
00044   r *= r;
00045   double r_weib = mean_sigma_sq(k);
00046   double res = r_weib-r;
00047   gradient[0]=2*res*neu/den;
00048 }
00049 
00050 double bsta_weibull_cost_function::lambda(double k) const
00051 {
00052   assert(k>0);
00053   return mean_/vnl_gamma((1+k)/k);
00054 }
00055 #endif // bsta_fit_weibull_cxx_