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 
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 
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 
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 
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 
186 
187 private:
191  GumbelMaxLikelihoodFitter & operator=(const GumbelMaxLikelihoodFitter & rhs);
192  };
193  }
194 }
195 
OpenMS::Math::GumbelMaxLikelihoodFitter
Implements a fitter for the Gumbel distribution.
Definition: GumbelMaxLikelihoodFitter.h:59
OpenMS::Math::GumbelMaxLikelihoodFitter::Functor::m_values
int m_values
Definition: GumbelMaxLikelihoodFitter.h:104
OpenMS::Math::GumbelMaxLikelihoodFitter::GumbelDistributionFitResult
struct to represent the parameters of a gumbel distribution
Definition: GumbelMaxLikelihoodFitter.h:65
double
DPosition.h
OpenMS::Math::GumbelMaxLikelihoodFitter::Functor
Definition: GumbelMaxLikelihoodFitter.h:93
OpenMS::Math::GumbelMaxLikelihoodFitter::GumbelDistributionFunctor::GumbelDistributionFunctor
GumbelDistributionFunctor(const std::vector< double > &data, const std::vector< double > &weights)
Definition: GumbelMaxLikelihoodFitter.h:117
OpenMS::Math::GumbelMaxLikelihoodFitter::Functor::InputType
Eigen::Matrix< Scalar, InputsAtCompileTime, 1 > InputType
Definition: GumbelMaxLikelihoodFitter.h:100
OpenMS::Math::GumbelMaxLikelihoodFitter::GumbelDistributionFunctor::operator()
int operator()(const Eigen::VectorXd &x, Eigen::VectorXd &fvec) const
Definition: GumbelMaxLikelihoodFitter.h:123
OpenMS::Math::GumbelMaxLikelihoodFitter::Functor::JacobianType
Eigen::Matrix< Scalar, ValuesAtCompileTime, InputsAtCompileTime > JacobianType
Definition: GumbelMaxLikelihoodFitter.h:102
OpenMS::Math::GumbelMaxLikelihoodFitter::GumbelDistributionFunctor::m_weights
const std::vector< double > & m_weights
Definition: GumbelMaxLikelihoodFitter.h:140
OpenMS::Math::GumbelMaxLikelihoodFitter::Functor::inputs
int inputs() const
Definition: GumbelMaxLikelihoodFitter.h:109
OpenMS::Math::GumbelMaxLikelihoodFitter::Functor::Functor
Functor(int inputs, int values)
Definition: GumbelMaxLikelihoodFitter.h:107
OpenMS
Main OpenMS namespace.
Definition: FeatureDeconvolution.h:46
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::Functor::ValueType
Eigen::Matrix< Scalar, ValuesAtCompileTime, 1 > ValueType
Definition: GumbelMaxLikelihoodFitter.h:101
OpenMS::Math::GumbelMaxLikelihoodFitter::GumbelDistributionFitResult::b
double b
scale parameter b
Definition: GumbelMaxLikelihoodFitter.h:76
OpenMS::Math::GumbelMaxLikelihoodFitter::GumbelDistributionFunctor::m_data
const std::vector< double > & m_data
Definition: GumbelMaxLikelihoodFitter.h:139
OpenMS::Math::GumbelMaxLikelihoodFitter::GumbelDistributionFitResult::GumbelDistributionFitResult
GumbelDistributionFitResult(double local_a, double local_b)
Definition: GumbelMaxLikelihoodFitter.h:67
OpenMS::Math::GumbelMaxLikelihoodFitter::GumbelDistributionFunctor
Definition: GumbelMaxLikelihoodFitter.h:114
OpenMS::Math::GumbelMaxLikelihoodFitter::Functor::values
int values() const
Definition: GumbelMaxLikelihoodFitter.h:110
OpenMS::Math::GumbelMaxLikelihoodFitter::GumbelDistributionFitResult::a
double a
location parameter a
Definition: GumbelMaxLikelihoodFitter.h:74
OpenMS::Math::GumbelMaxLikelihoodFitter::Functor::Functor
Functor()
Definition: GumbelMaxLikelihoodFitter.h:106
String.h
OpenMS::Math::GumbelMaxLikelihoodFitter::init_param_
GumbelDistributionFitResult init_param_
Definition: GumbelMaxLikelihoodFitter.h:185
OpenMS::Math::GumbelMaxLikelihoodFitter::Functor::Scalar
_Scalar Scalar
Definition: GumbelMaxLikelihoodFitter.h:95
OpenMS::Exception::UnableToFit
Exception used if an error occurred while fitting a model to a given dataset.
Definition: Exception.h:676