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-2017.
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: Timo Sachsenberg $
32 // $Authors: Christian Ehrlich, Chris Bielow $
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  /*
53  @brief Estimates model parameters for a quadratic equation
54 
55  The quadratic equation is of the form
56  y = a + b*x + c*x^2
57 
58  Weighted inputs are supported.
59 
60  */
61  class OPENMS_DLLAPI QuadraticRegression
62  {
63 public:
65 
67  template <typename Iterator>
68  void computeRegression(Iterator x_begin, Iterator x_end, Iterator y_begin);
69 
71  template <typename Iterator>
72  void computeRegressionWeighted(Iterator x_begin, Iterator x_end, Iterator y_begin, Iterator w_begin);
73 
75  double eval(double x) const;
76 
78  static double eval(double A, double B, double C, double x);
79 
80  double getA() const;
81  double getB() const;
82  double getC() const;
83  double getChiSquared() const;
84 
85 protected:
86  double a_;
87  double b_;
88  double c_;
89  double chi_squared_;
90  }; //class
91 
92  namespace
93  {
94  //x, y must be of same size
95  template <typename Iterator>
96  double computeChiSquareWeighted_(
97  Iterator x_begin, const Iterator& x_end, Iterator y_begin, Iterator w_begin,
98  const double a, const double b, const double c)
99  {
100  double chi_squared(0.0);
101  for (; x_begin != x_end; ++x_begin, ++y_begin, ++w_begin)
102  {
103  const double& x = *x_begin;
104  const double& y = *y_begin;
105  const double& weight = *w_begin;
106  chi_squared += weight * std::pow(y - a - b * x - c * x * x, 2);
107  }
108 
109  return chi_squared;
110  }
111 
112  } // anonymous namespace
113 
114  template <typename Iterator>
115  void QuadraticRegression::computeRegression(Iterator x_begin, Iterator x_end, Iterator y_begin)
116  {
117  std::vector<double> weights(std::distance(x_begin, x_end), 1);
118  computeRegressionWeighted<Iterator>(x_begin, x_end, y_begin, weights.begin());
119  }
120 
121  template <typename Iterator>
123  Iterator x_begin, Iterator x_end, Iterator y_begin, Iterator w_begin)
124  {
125  // Compute the linear fit of a quadratic function.
126  // Get the coefficients for y = w_1*a +w_2*b*x + w_3*c*x^2.
127  std::vector<Wm5::Vector2d> points = iteratorRange2Wm5Vectors(x_begin, x_end, y_begin);
128  // Compute sums for linear system. copy&paste from GeometricTools Wm5ApprLineFit2.cpp
129  // and modified to allow quadratic functions
130  int numPoints = static_cast<Int>(points.size());
131  double sumX = 0, sumXX = 0, sumXXX = 0, sumXXXX = 0;
132  double sumY = 0, sumXY = 0, sumXXY = 0;
133  double sumW = 0;
134 
135  Iterator wIter = w_begin;
136  for (int i = 0; i < numPoints; ++i, ++wIter)
137  {
138 
139  double x = points[i].X();
140  double y = points[i].Y();
141  double weight = *wIter;
142 
143  sumX += weight * x;
144  sumXX += weight * x * x;
145  sumXXX += weight * x * x * x;
146  sumXXXX += weight * x * x * x * x;
147 
148  sumY += weight * y;
149  sumXY += weight * x * y;
150  sumXXY += weight * x * x * y;
151 
152  sumW += weight;
153  }
154  //create matrices to solve Ax = B
155  double A[3][3] =
156  {
157  {sumW, sumX, sumXX},
158  {sumX, sumXX, sumXXX},
159  {sumXX, sumXXX, sumXXXX}
160  };
161  double B[3] =
162  {
163  sumY,
164  sumXY,
165  sumXXY
166  };
167  double X[3] = {0, 0, 0};
168 
169  bool nonsingular = Wm5::LinearSystem<double>().Solve3(A, B, X);
170  if (nonsingular)
171  {
172  a_ = X[0];
173  b_ = X[1];
174  c_ = X[2];
175  chi_squared_ = computeChiSquareWeighted_(x_begin, x_end, y_begin, w_begin, a_, b_, c_);
176  }
177  else
178  {
179  throw Exception::UnableToFit(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "UnableToFit-QuadraticRegression", "Could not fit a linear model to the data");
180  }
181  }
182 
183  } //namespace
184 } //namespace
185 
186 #endif // OPENMS_MATH_STATISTICS_QUADRATICREGRESSION_H
double c_
Definition: QuadraticRegression.h:88
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:61
double chi_squared_
Definition: QuadraticRegression.h:89
double b_
Definition: QuadraticRegression.h:87
void computeRegression(Iterator x_begin, Iterator x_end, Iterator y_begin)
Definition: QuadraticRegression.h:115
double a_
Definition: QuadraticRegression.h:86
void computeRegressionWeighted(Iterator x_begin, Iterator x_end, Iterator y_begin, Iterator w_begin)
Definition: QuadraticRegression.h:122
Exception used if an error occurred while fitting a model to a given dataset.
Definition: Exception.h:677
int Int
Signed integer type.
Definition: Types.h:103

OpenMS / TOPP release 2.3.0 Documentation generated on Tue Jan 9 2018 18:22:03 using doxygen 1.8.13