mxlib
c++ tools for analyzing astronomical data and other tasks by Jared R. Males. [git repo]
fitGammaDistribution.hpp
Go to the documentation of this file.
1 /** \file fitGammaDistribution.hpp
2  * \author Jared R. Males
3  * \brief Tools for fitting the GammaDistribution distribution to data.
4  * \ingroup fitting_files
5  *
6  */
7 
8 //***********************************************************************//
9 // Copyright 2023 Jared R. Males (jaredmales@gmail.com)
10 //
11 // This file is part of mxlib.
12 //
13 // mxlib is free software: you can redistribute it and/or modify
14 // it under the terms of the GNU General Public License as published by
15 // the Free Software Foundation, either version 3 of the License, or
16 // (at your option) any later version.
17 //
18 // mxlib is distributed in the hope that it will be useful,
19 // but WITHOUT ANY WARRANTY; without even the implied warranty of
20 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 // GNU General Public License for more details.
22 //
23 // You should have received a copy of the GNU General Public License
24 // along with mxlib. If not, see <http://www.gnu.org/licenses/>.
25 //***********************************************************************//
26 
27 #ifndef fitGammaDistribution_hpp
28 #define fitGammaDistribution_hpp
29 
30 #include "levmarInterface.hpp"
31 #include "../func/gammaDistribution.hpp"
32 
33 namespace mx
34 {
35 namespace math
36 {
37 namespace fit
38 {
39 
40 template<typename realT>
41 struct array2FitGammaDistribution;
42 
43 
44 /** \defgroup gammaDist_peak_fit Gamma Distribution
45  * \brief Fitting the Gamma Distribution to data.
46  *
47  * The Gamma Distribution is fit to data.
48  *
49  * \ingroup peak_fit
50  */
51 
52 ///Class to manage fitting the GammaDistribution Distribution to data via the \ref levmarInterface
53 /** In addition to the requirements on fitterT specified by \ref levmarInterface
54  * this class also requires this definition in fitterT
55  * \code
56  * static const int nparams = 3;
57  * \endcode
58  * where the number 3 is replaced by the number of parameters that fitterT expects to fit.
59  *
60  * \tparam fitterT a type meeting the above requirements.
61  *
62  * \ingroup gammaDist_peak_fit
63  *
64  */
65 template<typename fitterT>
66 class fitGammaDistribution : public levmarInterface<fitterT>
67 {
68 
69 public:
70 
71  typedef typename fitterT::realT realT;
72 
73  static const int nparams = fitterT::nparams;
74 
76 
77  void initialize()
78  {
79  this->allocate_params(nparams);
80  this->adata = &arr;
81  }
82 
84  {
85  initialize();
86  }
87 
89  {
90  }
91 
92  ///Set the initial guess.
93  void setGuess( realT x0, ///< [in] the location parameter
94  realT k, ///< [in] the shape parameter
95  realT theta, ///< [in] the scale parameter
96  realT denom ///< [in] the denominator or 1/peak-scale
97  )
98  {
99  static_assert( nparams==4 , "fitGammaDistribution: Wrong setGuess called for no location parameter.");
100 
101  this->p[2] = x0; //this is p[2] to make it easy to leave out
102  this->p[0] = k;
103  this->p[1] = theta;
104  this->p[3] = denom;
105  }
106 
107  ///Set the initial guess.
108  void setGuess( realT x0, ///< [in] the location parameter
109  realT k, ///< [in] the shape parameter
110  realT theta ///< [in] the scale parameter
111  )
112  {
113  static_assert( nparams==3 , "fitGammaDistribution: Wrong setGuess called for no location parameter.");
114 
115  this->p[2] = x0; //this is p[2] to make it easy to leave out
116  this->p[0] = k;
117  this->p[1] = theta;
118  }
119 
120  ///Set the initial guess when no location parameter is used.
121  void setGuess( realT k, ///< [in] the shape parameter
122  realT theta ///< [in] the scale parameter
123  )
124  {
125  static_assert( nparams==2 , "fitGammaDistribution: Wrong setGuess called for location parameter.");
126 
127  this->p[0] = k;
128  this->p[1] = theta;
129  }
130 
131  void setArray(realT *data, int n)
132  {
133  arr.data = data;
134  arr.n = n;
135 
136  this->n = n;
137 
138  }
139 
140  void x0( realT nx0 )
141  {
142  arr.x0 = nx0;
143  if(nparams == 3)
144  {
145  this->p[2] = nx0; //this is p[2] to make it easy to leave out
146  }
147  }
148 
149  void k( realT nk )
150  {
151  arr.k = nk;
152  this->p[0] = nk;
153  }
154 
155  void theta( realT na )
156  {
157  arr.lambda = na;
158  this->p[1] = na;
159  }
160 
161  void denom( realT nd)
162  {
163  arr.denom = nd;
164  if(nparams == 4)
165  {
166  this->p[3] = nd;
167  }
168  }
169 
170  int fit()
171  {
173  }
174 
175  realT x0()
176  {
177  if(nparams == 3)
178  {
179  return this->p[2];
180  }
181  else
182  {
183  return 0;
184  }
185  }
186 
187  realT k()
188  {
189  return this->p[0];
190  }
191 
192  realT theta()
193  {
194  return this->p[1];
195  }
196 
197  realT denom()
198  {
199  if(nparams == 4)
200  {
201  return this->p[3];
202  }
203  else
204  {
205  return func::gammaDistributionDenom<realT>(this->p[0], this->p[1]);
206  }
207  }
208 
209 
210 };
211 
212 
213 ///Wrapper for a native array to pass to \ref levmarInterface, with GammaDistribution details.
214 /** \ingroup gammaDist_peak_fit
215  */
216 template<typename realT>
218 {
219  realT * data {nullptr}; ///< Pointer to the array
220  size_t n {0}; ///< dimension of the array
221 
222  realT x0 {0}; ///< the location parameter.
223  realT k {0}; ///< the shape parameter
224  realT theta {0}; ///< the scale parameter
225  realT denom {0}; ///< the denominator or 1/peak-scale
226 };
227 
228 ///\ref levmarInterface fitter structure for the shifted Gamma Distribution with arbitrary peak scaling
229 /**
230  *
231  * \ingroup gammaDist_peak_fit
232  *
233  */
234 template<typename _realT>
236 {
237  typedef _realT realT;
238 
239  static const int nparams = 4;
240 
241  static void func(realT *p, realT *hx, int m, int n, void *adata)
242  {
244 
245  for(int i=0;i<arr->n; i++)
246  {
247  hx[i] = func::gammaDistribution<realT>(i, p[2], p[0], p[1], p[3]) - arr->data[i];
248  }
249  }
250 };
251 
252 ///\ref levmarInterface fitter structure for the shifted Gamma Distribution
253 /**
254  *
255  * \ingroup gammaDist_peak_fit
256  *
257  */
258 template<typename _realT>
260 {
261  typedef _realT realT;
262 
263  static const int nparams = 3;
264 
265  static void func(realT *p, realT *hx, int m, int n, void *adata)
266  {
268 
269  realT denom = func::gammaDistributionDenom<realT>(p[0], p[1]);
270 
271  for(int i=0;i<arr->n; i++)
272  {
273  hx[i] = func::gammaDistribution<realT>(i, p[2], p[0], p[1], denom) - arr->data[i];
274  }
275  }
276 };
277 
278 ///\ref levmarInterface fitter structure for the Gamma Distribution
279 /**
280  *
281  * \ingroup gammaDist_peak_fit
282  *
283  */
284 template<typename _realT>
286 {
287  typedef _realT realT;
288 
289  static const int nparams = 2;
290 
291  static void func(realT *p, realT *hx, int m, int n, void *adata)
292  {
294 
295  realT denom = func::gammaDistributionDenom<realT>(p[0], p[1]);
296 
297  for(int i=0;i<arr->n; i++)
298  {
299  hx[i] = func::gammaDistribution<realT>(i, 0.0, p[0], p[1], denom) - arr->data[i];
300  }
301  }
302 };
303 
304 
305 } //namespace fit
306 } //namespace math
307 } //namespace mx
308 
309 #endif //fitGammaDistribution_hpp
310 
Class to manage fitting the GammaDistribution Distribution to data via the levmarInterface.
void setGuess(realT x0, realT k, realT theta, realT denom)
Set the initial guess.
void setGuess(realT x0, realT k, realT theta)
Set the initial guess.
void setGuess(realT k, realT theta)
Set the initial guess when no location parameter is used.
A templatized interface to the levmar package.
void allocate_params()
Allocate parameters array based on previous call to nParams.
realT * p
Parameter array. On input is the initial estimates. On output has the estimated solution.
int n
I: measurement vector dimension.
void * adata
Pointer to possibly additional data, passed uninterpreted to func & jacf.
A c++ interface to the templatized levmar minimization routines..
The mxlib c++ namespace.
Definition: mxError.hpp:107
Wrapper for a native array to pass to levmarInterface, with GammaDistribution details.
realT denom
the denominator or 1/peak-scale
levmarInterface fitter structure for the Gamma Distribution
levmarInterface fitter structure for the shifted Gamma Distribution
levmarInterface fitter structure for the shifted Gamma Distribution with arbitrary peak scaling