← Back to branch summary

~ubuntu-branches/ubuntu/wily/openms/wily

« back to all changes in this revision

Viewing changes to source/MATH/STATISTICS/GumbelDistributionFitter.C

  • Committer: Package Import Robot
  • Author(s): Filippo Rusconi
  • Date: 2012-11-12 15:58:12 UTC
  • Revision ID: package-import@ubuntu.com-20121112155812-vr15wtg9b50cuesg
Tags: upstream-1.9.0
Import upstream version 1.9.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
source/MATH/STATISTICS/GumbelDistributionFitter.C
 
1
// -*- mode: C++; tab-width: 2; -*-
 
2
// vi: set ts=2:
 
3
//
 
4
// --------------------------------------------------------------------------
 
5
//                   OpenMS Mass Spectrometry Framework
 
6
// --------------------------------------------------------------------------
 
7
//  Copyright (C) 2003-2011 -- Oliver Kohlbacher, Knut Reinert
 
8
//
 
9
//  This library is free software; you can redistribute it and/or
 
10
//  modify it under the terms of the GNU Lesser General Public
 
11
//  License as published by the Free Software Foundation; either
 
12
//  version 2.1 of the License, or (at your option) any later version.
 
13
//
 
14
//  This library is distributed in the hope that it will be useful,
 
15
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
 
16
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 
17
//  Lesser General Public License for more details.
 
18
//
 
19
//  You should have received a copy of the GNU Lesser General Public
 
20
//  License along with this library; if not, write to the Free Software
 
21
//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
22
//
 
23
// --------------------------------------------------------------------------
 
24
// $Maintainer: David Wojnar $
 
25
// $Authors: David Wojnar $
 
26
// --------------------------------------------------------------------------
 
27
//
 
28
 
 
29
#include <sstream>
 
30
 
 
31
#include <cmath>
 
32
 
 
33
#include <gsl/gsl_sf_psi.h>
 
34
#include <OpenMS/MATH/STATISTICS/GumbelDistributionFitter.h>
 
35
 
 
36
using namespace std;
 
37
 
 
38
#define GUMBEL_DISTRIBUTION_FITTER_VERBOSE
 
39
#undef  GUMBEL_DISTRIBUTION_FITTER_VERBOSE
 
40
 
 
41
#ifdef GUMBEL_DISTRIBUTION_FITTER_VERBOSE
 
42
        #include <iostream>
 
43
        #include <gsl/gsl_rng.h>
 
44
        #include <gsl/gsl_randist.h>
 
45
        #include <gsl/gsl_vector.h>
 
46
        #include <gsl/gsl_blas.h>
 
47
        #include <gsl/gsl_multifit_nlin.h>
 
48
#endif
 
49
 
 
50
namespace OpenMS
 
51
{
 
52
        namespace Math
 
53
        {
 
54
                GumbelDistributionFitter::GumbelDistributionFitter()
 
55
                {
 
56
                        init_param_.a = 0.25;
 
57
                        init_param_.b = 0.1;
 
58
                }
 
59
        
 
60
                GumbelDistributionFitter::~GumbelDistributionFitter()
 
61
                {
 
62
                }
 
63
        
 
64
                void GumbelDistributionFitter::setInitialParameters(const GumbelDistributionFitResult& param)
 
65
                {
 
66
                        init_param_.a = param.a;
 
67
                        init_param_.b = param.b;
 
68
                }
 
69
        
 
70
                const String& GumbelDistributionFitter::getGnuplotFormula() const
 
71
                {
 
72
                        return gnuplot_formula_;
 
73
                }
 
74
        
 
75
                int GumbelDistributionFitter::gumbelDistributionFitterf_(const gsl_vector* x, void* params, gsl_vector* f)
 
76
                {
 
77
                        vector<DPosition<2> >* data = static_cast<vector<DPosition<2> >*>(params);
 
78
                        
 
79
                        double a = gsl_vector_get (x, 0);
 
80
                        double b = gsl_vector_get (x, 1);
 
81
                
 
82
                        UInt i = 0;
 
83
                        for (vector<DPosition<2> >::iterator it = data->begin(); it != data->end(); ++it)
 
84
                        {
 
85
                                double the_x = it->getX();
 
86
                                double z = exp((a - the_x)/b);
 
87
                                gsl_vector_set(f, i++, (z*exp(-1* z))/b - it->getY());
 
88
                        }
 
89
                
 
90
                        return GSL_SUCCESS;
 
91
                }
 
92
        
 
93
                // compute Jacobian matrix for the different parameters
 
94
                int GumbelDistributionFitter::gumbelDistributionFitterdf_(const gsl_vector* x, void* params, gsl_matrix* J)
 
95
                {
 
96
                        vector<DPosition<2> >* data = static_cast<vector<DPosition<2> >*>(params);
 
97
        
 
98
                        double a = gsl_vector_get (x, 0);
 
99
                        double b = gsl_vector_get (x, 1);
 
100
                        UInt i(0);
 
101
                        for (vector<DPosition<2> >::iterator it = data->begin(); it != data->end(); ++it, ++i)
 
102
                        {
 
103
                                double the_x = it->getX();
 
104
                                double z = exp((a - the_x)/b);
 
105
                                double f = z* exp (-1 * z);
 
106
                                double part_dev_a = (f - pow(z,2) * exp(-1*z))/pow(b,2);
 
107
                                gsl_matrix_set(J, i, 0,part_dev_a);
 
108
                                double dev_z =  ((the_x -a)/pow(b,2));
 
109
                                double cum = f*dev_z;
 
110
                                double part_dev_b = (( cum - z* cum)* b - f)/pow(b,2);
 
111
                                gsl_matrix_set(J,i,1,part_dev_b);
 
112
                                
 
113
 
 
114
                                
 
115
 
 
116
                        }
 
117
                  return GSL_SUCCESS;
 
118
                }
 
119
        
 
120
                int GumbelDistributionFitter::gumbelDistributionFitterfdf_(const gsl_vector* x, void* params, gsl_vector* f, gsl_matrix* J)
 
121
                {
 
122
                  gumbelDistributionFitterf_(x, params, f);
 
123
                  gumbelDistributionFitterdf_(x, params, J);
 
124
                  return GSL_SUCCESS;
 
125
                }
 
126
        
 
127
#ifdef GUMBEL_DISTRIBUTION_FITTER_VERBOSE
 
128
                void GumbelDistributionFitter::printState_(size_t iter, gsl_multifit_fdfsolver * s)
 
129
                {
 
130
                  printf ("iter: %3u x = % 15.8f % 15.8f "
 
131
                          "|f(x)| = %g\n",
 
132
                          (unsigned int)iter,
 
133
                          gsl_vector_get(s->x, 0), 
 
134
                          gsl_vector_get(s->x, 1),
 
135
                          gsl_blas_dnrm2(s->f));
 
136
                }
 
137
#endif
 
138
        
 
139
                GumbelDistributionFitter::GumbelDistributionFitResult GumbelDistributionFitter::fit(vector<DPosition<2> >& input)
 
140
                {
 
141
                  const gsl_multifit_fdfsolver_type* T = NULL;
 
142
                  gsl_multifit_fdfsolver* s = NULL;
 
143
                
 
144
                  int status = 0;
 
145
                  size_t iter = 0;
 
146
                
 
147
                  const size_t p = 2;
 
148
                
 
149
                  gsl_multifit_function_fdf f;
 
150
                  double x_init[2] = { init_param_.a, init_param_.b };
 
151
                  gsl_vector_view x = gsl_vector_view_array (x_init, p);
 
152
                  const gsl_rng_type * type = NULL;
 
153
                  gsl_rng* r = NULL;
 
154
                
 
155
                  gsl_rng_env_setup();
 
156
                
 
157
                  type = gsl_rng_default;
 
158
                  r = gsl_rng_alloc (type);
 
159
                
 
160
                  f.f = &gumbelDistributionFitterf_;
 
161
                  f.df = &gumbelDistributionFitterdf_;
 
162
                  f.fdf = &gumbelDistributionFitterfdf_;
 
163
                  f.n = input.size();
 
164
                  f.p = p;
 
165
                  f.params = &input;
 
166
                
 
167
                  T = gsl_multifit_fdfsolver_lmsder;
 
168
                  s = gsl_multifit_fdfsolver_alloc (T, input.size(), p);
 
169
                  gsl_multifit_fdfsolver_set (s, &f, &x.vector);
 
170
        
 
171
                        #ifdef GUMBEL_DISTRIBUTION_FITTER_VERBOSE
 
172
                  printState_(iter, s);
 
173
                        #endif
 
174
                
 
175
                  do
 
176
                  {
 
177
                    ++iter;
 
178
                    status = gsl_multifit_fdfsolver_iterate (s);
 
179
                
 
180
                                #ifdef GUMBEL_DISTRIBUTION_FITTER_VERBOSE
 
181
                    printf ("status = %s\n", gsl_strerror (status));
 
182
                    printState_(iter, s);
 
183
                                #endif
 
184
                                if (status)
 
185
                                {
 
186
                            break;
 
187
                                }
 
188
        
 
189
                    status = gsl_multifit_test_delta (s->dx, s->x, 1e-4, 1e-4);
 
190
                                #ifdef GUMBEL_DISTRIBUTION_FITTER_VERBOSE
 
191
                                printf("Status = '%s'\n",  gsl_strerror(status));
 
192
                                #endif
 
193
                  }
 
194
                  while (status == GSL_CONTINUE && iter < 1000);
 
195
 
 
196
                        #ifdef GUMBEL_DISTRIBUTION_FITTER_VERBOSE
 
197
      printf("Final status = '%s'\n",  gsl_strerror(status));
 
198
      #endif
 
199
        
 
200
                        if (status!=GSL_SUCCESS)
 
201
                        {
 
202
                                gsl_rng_free(r);
 
203
                                gsl_multifit_fdfsolver_free(s);
 
204
        
 
205
                                throw Exception::UnableToFit(__FILE__,__LINE__,__PRETTY_FUNCTION__,"UnableToFit-GumbelDistributionFitter","Could not fit the gumbel distribution to the data");
 
206
                        }
 
207
                  
 
208
                        // write the result in a GumbelDistributionFitResult struct
 
209
                        GumbelDistributionFitResult result;
 
210
                        result.a = gsl_vector_get(s->x, 0);
 
211
                        result.b = gsl_vector_get(s->x, 1);
 
212
        
 
213
                        // build a formula with the fitted parameters for gnuplot
 
214
                        stringstream formula;
 
215
                        formula <<"f(x)=" << "(1/" << result.b <<") * " << "exp(( "<< result.a<< "- x)/"<< result.b <<") * exp(-exp(("<<result.a<<" - x)/"<<result.b<<"))";
 
216
                        gnuplot_formula_ = formula.str();
 
217
                        
 
218
#ifdef GUMBEL_DISTRIBUTION_FITTER_VERBOSE
 
219
                        cout << gnuplot_formula_ << endl;
 
220
#endif
 
221
                        
 
222
                        gsl_rng_free(r);
 
223
                        gsl_multifit_fdfsolver_free (s);
 
224
        
 
225
                        return result;
 
226
                }
 
227
 
 
228
        } //namespace Math
 
229
} // namespace OpenMS