GSL-- 0.11 - GSL--: fit/fit-ex.cc

00001 #include <gslmm/fit/fit.hh>
00002 #include <iostream>
00003 #include <sstream>
00004 #include <gslmm/random/gaussian.hh>
00005 #include <gslmm/error/error_handler.hh>
00006 #include <iomanip>
00007 #include <vector>
00008 
00009 struct expo 
00010 {
00011   std::vector<double> y;
00012   std::vector<double> sigma;
00013   
00014   int function(const gslmm::vector<double>& x, gslmm::vector<double>& f)
00015   {
00016     double A      = x[0];
00017     double lambda = x[1];
00018     double b      = x[2];
00019     for (size_t i = 0; i < y.size(); i++) {
00020       double t  = i;
00021       double yi = A * exp(-lambda * t) + b;
00022       f[i]      = (yi - y[i]) / sigma[i];
00023     }
00024     return GSL_SUCCESS;
00025   }
00026   int jacobian(const gslmm::vector<double>& x, gslmm::matrix<double>& j)
00027   {
00028     double A      = x[0];
00029     double lambda = x[1];
00030     double b      = x[2];
00031     for (size_t i = 0; i < y.size(); i++) {
00032       double t   = i;
00033       double s   = sigma[i];
00034       double e   = exp(-lambda * t);
00035       j(i, 0) = e / s;
00036       j(i, 1) = -t * A * e / s;
00037       j(i, 2) = 1 / s;
00038     }
00039     return GSL_SUCCESS;
00040   }
00041   int iterate() { return GSL_SUCCESS; }
00042   expo(size_t n) 
00043     : y(n), sigma(n)
00044   {
00045     gslmm::generator::env_setup();
00046     gslmm::generator g;
00047     gslmm::gaussian  r(0.1, g);
00048     
00049     for (size_t i = 0; i < n; i++) {
00050       double t = i;
00051       y[i]     = 1. + 5 * exp(-0.1 * t) + r.sample();
00052       sigma[i] = .1;
00053     }
00054   }
00055 };
00056 
00057 void 
00058 print_status(size_t iter, const gslmm::fitter<expo>& fit, bool cont) 
00059 {
00060   gslmm::vector<double> p(fit.position());
00061   gslmm::vector<double> f(fit.values());
00062   std::cout << std::setw(4) << iter << " x=("
00063             << std::setw(10) << p[0] << ","
00064             << std::setw(10) << p[1] << ","
00065             << std::setw(10) << p[2] << ") |f(x)|="
00066             << std::setw(10) << f.length() << "\t" 
00067             << (cont ? "continue" : "done") << std::endl;
00068 }
00069 
00070 int main()
00071 {
00072   expo e(40);
00073   gslmm::vector<double> initial(3);
00074   initial[0] = 1;
00075   gslmm::fitter<expo> fit;
00076   fit.initialize(e, gslmm::fitter<expo>::scaled_levenberg_marquardt, 
00077                initial, 40);
00078   
00079   print_status(0, fit, true);
00080   size_t iter = 1;
00081   int status  = 0;
00082   bool cont = true;
00083   do {
00084     status = fit.iterate();
00085     if (status) { 
00086       std::cout << "bad status " << status << std::endl;
00087       break;
00088     }
00089     cont = !fit.test_delta(1e-4, 1e-4);
00090     print_status(iter, fit, cont);
00091   } while (cont && ++iter < 500);
00092   gslmm::matrix<double> c(3,3);
00093   gslmm::vector<double> p(fit.position());
00094   fit.covariance(c, 0);
00095   
00096   for (size_t i = 0; i < 3; i++) {
00097     std::cout << (i == 0 ? "A" : (i == 1 ? "lambda" : "b"))
00098               << "\t\t= " << std::setw(10) << p[i] 
00099               << " +/" << std::setw(10) << sqrt(c(i,i)) << std::endl;
00100   }
00101   gslmm::vector<double> f(fit.values());
00102   double chisq = f.length();
00103   std::cout << "chi^2/dof\t= " << std::setw(10) << (chisq*chisq) 
00104             << " / 37 = " << (chisq*chisq)/37 << "\nstatus=" 
00105             << gslmm::error_handler<>::code_to_string(status) << std::endl;
00106   return status;
00107 }
00108 
00109   
00110 
00111   
00112