00001 #include <map>
00002 #include <stdexcept>
00003 #include <string>
00004 #include <math.h>
00005 template<typename Number>
00006 class LinearFit
00007 {
00008 public:
00009 typedef map<Number, Number> Points;
00010 LinearFit();
00011 LinearFit(Points * points);
00012 LinearFit(const LinearFit & fit);
00013 virtual ~LinearFit();
00014 const LinearFit & operator=(const LinearFit & fit);
00015 virtual void fit();
00016 virtual void fit(Points * points);
00017 const Number getSlope() const;
00018 const Number getIntercept() const;
00019 const Number getRegressionCoefficient() const;
00020 protected:
00021 Points * _points;
00023 Number _a;
00025 Number _b;
00027 Number _r;
00028 };
00029
00030 template<typename Number>
00031 LinearFit<Number>::LinearFit()
00032 : _points(0), _a(0), _b(0), _r(0)
00033 {}
00034
00035 template<typename Number>
00036 LinearFit<Number>::LinearFit(Points * points)
00037 : _points(points), _a(0), _b(0), _r(0)
00038 {}
00039
00040 template<typename Number>
00041 LinearFit<Number>::LinearFit(const LinearFit<Number> & fit)
00042 : _points(fit._points), _a(fit._a), _b(fit._b), _r(fit._r)
00043 {}
00044
00045 template<typename Number>
00046 const LinearFit<Number> & LinearFit<Number>::operator=(const LinearFit<Number> & fit)
00047 {
00048 _points = fit._points;
00049 _a = fit._a;
00050 _b = fit._b;
00051 _r = fit._r;
00052
00053 return *this;
00054 }
00055
00059 template<typename Number>
00060 void LinearFit<Number>::fit(Points * points)
00061 {
00062 _points = points;
00063 _a = _b = _r = 0;
00064 fit();
00065 }
00066
00067 template<typename Number>
00068 void LinearFit<Number>::fit()
00069 {
00070 int n = _points->size();
00071 if (n < 2)
00072 throw runtime_error("LinearFit<Number>::fit() - FATAL - n<2");
00073 Number sumx=0,sumy=0,sumx2=0,sumy2=0,sumxy=0;
00074 Number x, y,sxx,syy,sxy;
00075
00076 typename map<Number, Number>::const_iterator i;
00077 for (i = _points->begin(); i != _points->end(); i++)
00078 {
00079 x = i->first;
00080 y = i->second;
00081 sumx += x;
00082 sumy += y;
00083 sumx2 += (x*x);
00084 sumy2 += (y*y);
00085 sumxy += (x*y);
00086 }
00087 sxx = sumx2 - (sumx * sumx / n);
00088 syy = sumy2 - (sumy * sumy / n);
00089 sxy = sumxy - (sumx * sumy / n);
00090
00091 if (fabs(sxx) == 0)
00092 {
00093 _a = -99999.;
00094 _b = -99999.;
00095 _r = -99999.;
00096 return;
00097 }
00098
00099 _a = sxy/sxx;
00100 _b = sumy/n - _a*sumx/n;
00101
00102 if (fabs(syy) == 0)
00103 _r = 1;
00104 else
00105 _r = sxy/::sqrt(sxx*syy);
00106 }
00107
00108 template<typename Number>
00109 LinearFit<Number>::~LinearFit()
00110 {}
00111
00112 template<typename Number>
00113 inline const Number LinearFit<Number>::getSlope() const
00114 {
00115 return _a;
00116 }
00117
00118 template<typename Number>
00119 inline const Number LinearFit<Number>::getIntercept() const
00120 {
00121 return _b;
00122 }
00123
00124 template<typename Number>
00125 inline const Number LinearFit<Number>::getRegressionCoefficient() const
00126 {
00127 return _r;
00128 }
00129
