| OpenMS 2.5.0

GumbelMaxLikelihoodFitter.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: Julianus Pfeuffer $
32 // $Authors: Julianus Pfeuffer $
33 // --------------------------------------------------------------------------
34 //
35 #pragma once
36
37 #include <OpenMS/DATASTRUCTURES/String.h>
38 #include <OpenMS/DATASTRUCTURES/DPosition.h>
39 #include <unsupported/Eigen/NonLinearOptimization>
40 #include <vector>
41
42
43 namespace OpenMS
44 {
45 namespace Math
46 {
59 class OPENMS_DLLAPI GumbelMaxLikelihoodFitter
60 {
61
62 public:
63
65 struct GumbelDistributionFitResult
66 {
67 GumbelDistributionFitResult(double local_a, double local_b) :
68 a(local_a),
69 b(local_b)
70 {
71 }
72
74 double a;
76 double b;
77
78 double log_eval_no_normalize(double x) const;
79 };
80
82 GumbelMaxLikelihoodFitter();
84 GumbelMaxLikelihoodFitter(GumbelDistributionFitResult init);
86 virtual ~GumbelMaxLikelihoodFitter();
87
89 void setInitialParameters(const GumbelDistributionFitResult & result);
90
91 // Generic functor
92 template<typename _Scalar, int NX = Eigen::Dynamic, int NY = Eigen::Dynamic>
93 struct Functor
94 {
95 typedef _Scalar Scalar;
96 enum {
97 InputsAtCompileTime = NX,
98 ValuesAtCompileTime = NY
99 };
100 typedef Eigen::Matrix<Scalar,InputsAtCompileTime,1> InputType;
101 typedef Eigen::Matrix<Scalar,ValuesAtCompileTime,1> ValueType;
102 typedef Eigen::Matrix<Scalar,ValuesAtCompileTime,InputsAtCompileTime> JacobianType;
103
104 int m_inputs, m_values;
105
106 Functor() : m_inputs(InputsAtCompileTime), m_values(ValuesAtCompileTime) {}
107 Functor(int inputs, int values) : m_inputs(inputs), m_values(values) {}
108
109 int inputs() const { return m_inputs; }
110 int values() const { return m_values; }
111
112 };
113
114 struct GumbelDistributionFunctor : Functor<double>
115 {
116
117 GumbelDistributionFunctor(const std::vector<double>& data, const std::vector<double>& weights):
118 Functor<double>(2,2),
119 m_data(data), m_weights(weights)
120 {
121 }
122
123 int operator()(const Eigen::VectorXd &x, Eigen::VectorXd &fvec) const
124 {
125 fvec(0) = 0.0;
126 double sigma = fabs(x(1));
127 double logsigma = log(sigma);
128 auto wit = m_weights.cbegin();
129 for (auto it = m_data.cbegin(); it != m_data.cend(); ++it, ++wit)
130 {
131 double diff = (*it - x(0)) / sigma;
132 fvec(0) += *wit * (-logsigma - diff - exp(-diff));
133 }
134 double foo = -fvec(0);
135 fvec(0) = foo;
136 fvec(1) = 0.0;
137 return 0;
138 }
139 const std::vector<double>& m_data;
140 const std::vector<double>& m_weights;
141 };
142
152 GumbelDistributionFitResult fitWeighted(const std::vector<double> & x, const std::vector<double> & w)
153 {
154
155 Eigen::VectorXd x_init (2);
156 x_init(0) = init_param_.a;
157 x_init(1) = init_param_.b;
158 GumbelDistributionFunctor functor (x, w);
159 Eigen::NumericalDiff<GumbelDistributionFunctor> numDiff(functor);
160 Eigen::LevenbergMarquardt<Eigen::NumericalDiff<GumbelDistributionFunctor>,double> lm(numDiff);
161 Eigen::LevenbergMarquardtSpace::Status status = lm.minimize(x_init);
162
163 //the states are poorly documented. after checking the source, we believe that
164 //all states except NotStarted, Running and ImproperInputParameters are good
165 //termination states.
166 if (status <= Eigen::LevenbergMarquardtSpace::Status::ImproperInputParameters)
167 {
168 throw Exception::UnableToFit(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "UnableToFit-GumbelMaxLikelihoodFitter", "Could not fit the gumbel distribution to the data");
169 }
170
171 #ifdef GUMBEL_DISTRIBUTION_FITTER_VERBOSE
172 // build a formula with the fitted parameters for gnuplot
173 stringstream formula;
174 formula << "f(x)=" << "(1/" << x_init(1) << ") * " << "exp(( " << x_init(0) << "- x)/" << x_init(1) << ") * exp(-exp((" << x_init(0) << " - x)/" << x_init(1) << "))";
175 cout << formula.str() << endl;
176 #endif
177 init_param_.a = x_init(0);
178 init_param_.b = fabs(x_init(1));
179
180 return {x_init(0), fabs(x_init(1))};
181 }
182
183 protected:
184
185 GumbelDistributionFitResult init_param_;
186
187 private:
189 GumbelMaxLikelihoodFitter(const GumbelMaxLikelihoodFitter & rhs);
191 GumbelMaxLikelihoodFitter & operator=(const GumbelMaxLikelihoodFitter & rhs);
192 };
193 }
194 }
195
OpenMS::Math::GumbelMaxLikelihoodFitter::GumbelDistributionFitResult::b
double b
scale parameter b
Definition: GumbelMaxLikelihoodFitter.h:76
OpenMS::Math::GumbelMaxLikelihoodFitter::Functor::JacobianType
Eigen::Matrix< Scalar, ValuesAtCompileTime, InputsAtCompileTime > JacobianType
Definition: GumbelMaxLikelihoodFitter.h:102
OpenMS::Math::GumbelMaxLikelihoodFitter::Functor::m_values
int m_values
Definition: GumbelMaxLikelihoodFitter.h:104
OpenMS::Math::GumbelMaxLikelihoodFitter::GumbelDistributionFunctor::operator()
int operator()(const Eigen::VectorXd &x, Eigen::VectorXd &fvec) const
Definition: GumbelMaxLikelihoodFitter.h:123
OpenMS::Math::GumbelMaxLikelihoodFitter::GumbelDistributionFunctor::m_data
const std::vector< double > & m_data
Definition: GumbelMaxLikelihoodFitter.h:139
OpenMS::Math::GumbelMaxLikelihoodFitter::Functor::inputs
int inputs() const
Definition: GumbelMaxLikelihoodFitter.h:109
OpenMS::Math::GumbelMaxLikelihoodFitter::GumbelDistributionFitResult
struct to represent the parameters of a gumbel distribution
Definition: GumbelMaxLikelihoodFitter.h:65
OpenMS::Math::GumbelMaxLikelihoodFitter::Functor::Functor
Functor()
Definition: GumbelMaxLikelihoodFitter.h:106
OpenMS::Math::GumbelMaxLikelihoodFitter::GumbelDistributionFitResult::a
double a
location parameter a
Definition: GumbelMaxLikelihoodFitter.h:74
OpenMS::Math::GumbelMaxLikelihoodFitter::Functor::Functor
Functor(int inputs, int values)
Definition: GumbelMaxLikelihoodFitter.h:107
OpenMS::Exception::UnableToFit
Exception used if an error occurred while fitting a model to a given dataset.
Definition: Exception.h:676
OpenMS::Math::GumbelMaxLikelihoodFitter::GumbelDistributionFunctor
Definition: GumbelMaxLikelihoodFitter.h:114
OpenMS::Math::GumbelMaxLikelihoodFitter
Implements a fitter for the Gumbel distribution.
Definition: GumbelMaxLikelihoodFitter.h:59
String.h
DPosition.h
OpenMS::Math::GumbelMaxLikelihoodFitter::Functor::ValueType
Eigen::Matrix< Scalar, ValuesAtCompileTime, 1 > ValueType
Definition: GumbelMaxLikelihoodFitter.h:101
OpenMS
Main OpenMS namespace.
Definition: FeatureDeconvolution.h:46
OpenMS::Math::GumbelMaxLikelihoodFitter::Functor
Definition: GumbelMaxLikelihoodFitter.h:93
OpenMS::Math::GumbelMaxLikelihoodFitter::fitWeighted
GumbelDistributionFitResult fitWeighted(const std::vector< double > &x, const std::vector< double > &w)
Fits a gumbel distribution to the given data x values. Fills a weighted histogram first and generates...
Definition: GumbelMaxLikelihoodFitter.h:152
OpenMS::Math::GumbelMaxLikelihoodFitter::GumbelDistributionFunctor::m_weights
const std::vector< double > & m_weights
Definition: GumbelMaxLikelihoodFitter.h:140
Math
double
OpenMS::Math::GumbelMaxLikelihoodFitter::Functor::InputType
Eigen::Matrix< Scalar, InputsAtCompileTime, 1 > InputType
Definition: GumbelMaxLikelihoodFitter.h:100
OpenMS::Math::GumbelMaxLikelihoodFitter::Functor::values
int values() const
Definition: GumbelMaxLikelihoodFitter.h:110
OpenMS::Math::GumbelMaxLikelihoodFitter::GumbelDistributionFitResult::GumbelDistributionFitResult
GumbelDistributionFitResult(double local_a, double local_b)
Definition: GumbelMaxLikelihoodFitter.h:67
OpenMS::Math::GumbelMaxLikelihoodFitter::GumbelDistributionFunctor::GumbelDistributionFunctor
GumbelDistributionFunctor(const std::vector< double > &data, const std::vector< double > &weights)
Definition: GumbelMaxLikelihoodFitter.h:117
OpenMS::Math::GumbelMaxLikelihoodFitter::init_param_
GumbelDistributionFitResult init_param_
Definition: GumbelMaxLikelihoodFitter.h:185
OpenMS::Math::GumbelMaxLikelihoodFitter::Functor::Scalar
_Scalar Scalar
Definition: GumbelMaxLikelihoodFitter.h:95