Home  · Classes  · Annotated Classes  · Modules  · Members  · Namespaces  · Related Pages
QuadraticRegression.h
Go to the documentation of this file.
1 // --------------------------------------------------------------------------
2 // OpenMS -- Open-Source Mass Spectrometry
3 // --------------------------------------------------------------------------
4 // Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
5 // ETH Zurich, and Freie Universitaet Berlin 2002-2015.
6 //
7 // This software is released under a three-clause BSD license:
8 // * Redistributions of source code must retain the above copyright
9 // notice, this list of conditions and the following disclaimer.
10 // * Redistributions in binary form must reproduce the above copyright
11 // notice, this list of conditions and the following disclaimer in the
12 // documentation and/or other materials provided with the distribution.
13 // * Neither the name of any author or any participating institution
14 // may be used to endorse or promote products derived from this software
15 // without specific prior written permission.
16 // For a full list of authors, refer to the file AUTHORS.
17 // --------------------------------------------------------------------------
18 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
21 // ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
22 // INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
23 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
24 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
25 // OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
26 // WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
27 // OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
28 // ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29 //
30 // --------------------------------------------------------------------------
31 // $Maintainer: Christian Ehrlich $
32 // $Authors: Christian Ehrlich $
33 // --------------------------------------------------------------------------
34 
35 #ifndef OPENMS_MATH_STATISTICS_QUADRATICREGRESSION_H
36 #define OPENMS_MATH_STATISTICS_QUADRATICREGRESSION_H
37 
38 #include <OpenMS/CONCEPT/Types.h>
41 
42 #include "Wm5Vector2.h"
43 #include "Wm5LinearSystem.h"
44 #include <iterator>
45 
46 using Wm5::LinearSystem;
47 
48 namespace OpenMS
49 {
50  namespace Math
51  {
52  class OPENMS_DLLAPI QuadraticRegression
53  {
54 public:
56 
58  template <typename Iterator>
59  void computeRegression(Iterator x_begin, Iterator x_end,
60  Iterator y_begin);
61 
63  template <typename Iterator>
64  void computeRegressionWeighted(
65  Iterator x_begin, Iterator x_end, Iterator y_begin, Iterator w_begin);
66 
68  double eval(double x) const;
69 
70  double getA() const;
71  double getB() const;
72  double getC() const;
73  double getChiSquared() const;
74 
75 protected:
76  double a_;
77  double b_;
78  double c_;
79  double chi_squared_;
80  }; //class
81 
82  namespace
83  {
84  //x, y must be of same size
85  template <typename Iterator>
86  double computeChiSquareWeighted(
87  Iterator x_begin, Iterator x_end, Iterator y_begin, Iterator w_begin,
88  double a, double b, double c)
89  {
90  double chi_squared = 0.0;
91  Iterator xIter = x_begin;
92  Iterator yIter = y_begin;
93  Iterator wIter = w_begin;
94  for (; xIter != x_end; ++xIter, ++yIter, ++wIter)
95  {
96  double x = *xIter;
97  double y = *yIter;
98  double weigth = *wIter;
99  chi_squared += weigth * std::pow(y - a - b * x - c * x * x, 2);
100  }
101 
102  return chi_squared;
103  }
104 
105  }
106 
107  template <typename Iterator>
108  void QuadraticRegression::computeRegression(Iterator x_begin, Iterator x_end, Iterator y_begin)
109  {
110  std::vector<double> weights(std::distance(x_begin, x_end), 1);
111  computeRegressionWeighted<Iterator>(x_begin, x_end, y_begin, weights.begin());
112  }
113 
114  template <typename Iterator>
116  Iterator x_begin, Iterator x_end, Iterator y_begin, Iterator w_begin)
117  {
118  // Compute the linear fit of a quadratic function.
119  // Get the coefficients for y = w_1*a +w_2*b*x + w_3*c*x^2.
120  std::vector<Wm5::Vector2d> points = iteratorRange2Wm5Vectors(x_begin, x_end, y_begin);
121  // Compute sums for linear system. copy&paste from GeometricTools Wm5ApprLineFit2.cpp
122  // and modified to allow quadratic functions
123  int numPoints = static_cast<Int>(points.size());
124  double sumX = 0, sumXX = 0, sumXXX = 0, sumXXXX = 0;
125  double sumY = 0, sumXY = 0, sumXXY = 0;
126  double sumW = 0;
127 
128  Iterator wIter = w_begin;
129  for (int i = 0; i < numPoints; ++i, ++wIter)
130  {
131 
132  double x = points[i].X();
133  double y = points[i].Y();
134  double weight = *wIter;
135 
136  sumX += weight * x;
137  sumXX += weight * x * x;
138  sumXXX += weight * x * x * x;
139  sumXXXX += weight * x * x * x * x;
140 
141  sumY += weight * y;
142  sumXY += weight * x * y;
143  sumXXY += weight * x * x * y;
144 
145  sumW += weight;
146  }
147  //create matrices to solve Ax = B
148  double A[3][3] =
149  {
150  {sumW, sumX, sumXX},
151  {sumX, sumXX, sumXXX},
152  {sumXX, sumXXX, sumXXXX}
153  };
154  double B[3] =
155  {
156  sumY,
157  sumXY,
158  sumXXY
159  };
160  double X[3] = {0, 0, 0};
161 
162  bool nonsingular = Wm5::LinearSystem<double>().Solve3(A, B, X);
163  if (nonsingular)
164  {
165  a_ = X[0];
166  b_ = X[1];
167  c_ = X[2];
168  chi_squared_ = computeChiSquareWeighted(x_begin, x_end, y_begin, w_begin, a_, b_, c_);
169  }
170  else
171  {
172  throw Exception::UnableToFit(__FILE__, __LINE__, __PRETTY_FUNCTION__, "UnableToFit-QuadraticRegression", "Could not fit a linear model to the data");
173  }
174  }
175 
176  } //namespace
177 } //namespace
178 
179 #endif // OPENMS_MATH_STATISTICS_QUADRATICREGRESSION_H
double c_
Definition: QuadraticRegression.h:78
std::vector< Wm5::Vector2d > iteratorRange2Wm5Vectors(Iterator x_begin, Iterator x_end, Iterator y_begin)
Copies the distance(x_begin,x_end) elements starting at x_begin and y_begin into the Wm5::Vector...
Definition: RegressionUtils.h:45
const double c
Main OpenMS namespace.
Definition: FeatureDeconvolution.h:47
Definition: QuadraticRegression.h:52
double chi_squared_
Definition: QuadraticRegression.h:79
double b_
Definition: QuadraticRegression.h:77
void computeRegression(Iterator x_begin, Iterator x_end, Iterator y_begin)
Definition: QuadraticRegression.h:108
double a_
Definition: QuadraticRegression.h:76
void computeRegressionWeighted(Iterator x_begin, Iterator x_end, Iterator y_begin, Iterator w_begin)
Definition: QuadraticRegression.h:115
Exception used if an error occurred while fitting a model to a given dataset.
Definition: Exception.h:662
int Int
Signed integer type.
Definition: Types.h:96

OpenMS / TOPP release 2.0.0 Documentation generated on Fri May 29 2015 17:20:29 using doxygen 1.8.9.1