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