| OpenMS 2.5.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>
38 #include <OpenMS/CONCEPT/Exception.h>
39 #include <OpenMS/MATH/STATISTICS/RegressionUtils.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:
61 QuadraticRegression();
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>
112 void QuadraticRegression::computeRegression(Iterator x_begin, Iterator x_end, Iterator y_begin)
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>
119 void QuadraticRegression::computeRegressionWeighted(
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
OpenMS::Math::QuadraticRegression
Definition: QuadraticRegression.h:58
OpenMS::Exception::UnableToFit
Exception used if an error occurred while fitting a model to a given dataset.
Definition: Exception.h:676
Exception.h
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
Types.h
Iterator
OpenMS::Math::QuadraticRegression::computeRegressionWeighted
void computeRegressionWeighted(Iterator x_begin, Iterator x_end, Iterator y_begin, Iterator w_begin)
Definition: QuadraticRegression.h:119
OpenMS::Math::QuadraticRegression::a_
double a_
Definition: QuadraticRegression.h:83
OpenMS
Main OpenMS namespace.
Definition: FeatureDeconvolution.h:46
Math
OpenMS::Math::QuadraticRegression::chi_squared_
double chi_squared_
Definition: QuadraticRegression.h:86
OpenMS::Constants::c
const double c
int
OpenMS::Math::QuadraticRegression::b_
double b_
Definition: QuadraticRegression.h:84
OpenMS::Math::QuadraticRegression::computeRegression
void computeRegression(Iterator x_begin, Iterator x_end, Iterator y_begin)
Definition: QuadraticRegression.h:112
RegressionUtils.h
OpenMS::Math::QuadraticRegression::c_
double c_
Definition: QuadraticRegression.h:85