Line data Source code
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 : /** \defgroup gammaDist_peak_fit Gamma Distribution
44 : * \brief Fitting the Gamma Distribution to data.
45 : *
46 : * The Gamma Distribution is fit to data.
47 : *
48 : * \ingroup peak_fit
49 : */
50 :
51 : /// Class to manage fitting the GammaDistribution Distribution to data via the \ref levmarInterface
52 : /** In addition to the requirements on fitterT specified by \ref levmarInterface
53 : * this class also requires this definition in fitterT
54 : * \code
55 : * static const int nparams = 3;
56 : * \endcode
57 : * where the number 3 is replaced by the number of parameters that fitterT expects to fit.
58 : *
59 : * \tparam fitterT a type meeting the above requirements.
60 : *
61 : * \ingroup gammaDist_peak_fit
62 : *
63 : */
64 : template <typename fitterT>
65 : class fitGammaDistribution : public levmarInterface<fitterT>
66 : {
67 :
68 : public:
69 : typedef typename fitterT::realT realT;
70 :
71 : static const int nparams = fitterT::nparams;
72 :
73 : array2FitGammaDistribution<realT> arr;
74 :
75 : void initialize()
76 : {
77 : this->allocate_params( nparams );
78 : this->adata = &arr;
79 : }
80 :
81 : fitGammaDistribution()
82 : {
83 : initialize();
84 : }
85 :
86 : ~fitGammaDistribution()
87 : {
88 : }
89 :
90 : /// Set the initial guess.
91 : void setGuess( realT x0, ///< [in] the location parameter
92 : realT k, ///< [in] the shape parameter
93 : realT theta, ///< [in] the scale parameter
94 : realT denom ///< [in] the denominator or 1/peak-scale
95 : )
96 : {
97 : static_assert( nparams == 4, "fitGammaDistribution: Wrong setGuess called for no location parameter." );
98 :
99 : this->p[2] = x0; // this is p[2] to make it easy to leave out
100 : this->p[0] = k;
101 : this->p[1] = theta;
102 : this->p[3] = denom;
103 : }
104 :
105 : /// Set the initial guess.
106 : void setGuess( realT x0, ///< [in] the location parameter
107 : realT k, ///< [in] the shape parameter
108 : realT theta ///< [in] the scale parameter
109 : )
110 : {
111 : static_assert( nparams == 3, "fitGammaDistribution: Wrong setGuess called for no location parameter." );
112 :
113 : this->p[2] = x0; // this is p[2] to make it easy to leave out
114 : this->p[0] = k;
115 : this->p[1] = theta;
116 : }
117 :
118 : /// Set the initial guess when no location parameter is used.
119 : void setGuess( realT k, ///< [in] the shape parameter
120 : realT theta ///< [in] the scale parameter
121 : )
122 : {
123 : static_assert( nparams == 2, "fitGammaDistribution: Wrong setGuess called for location parameter." );
124 :
125 : this->p[0] = k;
126 : this->p[1] = theta;
127 : }
128 :
129 : void setArray( realT *data, int n )
130 : {
131 : arr.data = data;
132 : arr.n = n;
133 :
134 : this->n = n;
135 : }
136 :
137 : void x0( realT nx0 )
138 : {
139 : arr.x0 = nx0;
140 : if( nparams == 3 )
141 : {
142 : this->p[2] = nx0; // this is p[2] to make it easy to leave out
143 : }
144 : }
145 :
146 : void k( realT nk )
147 : {
148 : arr.k = nk;
149 : this->p[0] = nk;
150 : }
151 :
152 0 : void theta( realT na )
153 : {
154 0 : arr.theta = na;
155 0 : this->p[1] = na;
156 0 : }
157 :
158 : void denom( realT nd )
159 : {
160 : arr.denom = nd;
161 : if( nparams == 4 )
162 : {
163 : this->p[3] = nd;
164 : }
165 : }
166 :
167 : int fit()
168 : {
169 : return levmarInterface<fitterT>::fit();
170 : }
171 :
172 : realT x0()
173 : {
174 : if( nparams == 3 )
175 : {
176 : return this->p[2];
177 : }
178 : else
179 : {
180 : return 0;
181 : }
182 : }
183 :
184 : realT k()
185 : {
186 : return this->p[0];
187 : }
188 :
189 : realT theta()
190 : {
191 : return this->p[1];
192 : }
193 :
194 : realT denom()
195 : {
196 : if( nparams == 4 )
197 : {
198 : return this->p[3];
199 : }
200 : else
201 : {
202 : return func::gammaDistributionDenom<realT>( this->p[0], this->p[1] );
203 : }
204 : }
205 : };
206 :
207 : /// Wrapper for a native array to pass to \ref levmarInterface, with GammaDistribution details.
208 : /** \ingroup gammaDist_peak_fit
209 : */
210 : template <typename realT>
211 : struct array2FitGammaDistribution
212 : {
213 : realT *data{ nullptr }; ///< Pointer to the array
214 : size_t n{ 0 }; ///< dimension of the array
215 :
216 : realT x0{ 0 }; ///< the location parameter.
217 : realT k{ 0 }; ///< the shape parameter
218 : realT theta{ 0 }; ///< the scale parameter
219 : realT denom{ 0 }; ///< the denominator or 1/peak-scale
220 : };
221 :
222 : ///\ref levmarInterface fitter structure for the shifted Gamma Distribution with arbitrary peak scaling
223 : /**
224 : *
225 : * \ingroup gammaDist_peak_fit
226 : *
227 : */
228 : template <typename _realT>
229 : struct gammaDistribution_4param_fitter
230 : {
231 : typedef _realT realT;
232 :
233 : static const int nparams = 4;
234 :
235 : static void func( realT *p, realT *hx, int m, int n, void *adata )
236 : {
237 : array2FitGammaDistribution<realT> *arr = (array2FitGammaDistribution<realT> *)adata;
238 :
239 : for( int i = 0; i < arr->n; i++ )
240 : {
241 : hx[i] = func::gammaDistribution<realT>( i, p[2], p[0], p[1], p[3] ) - arr->data[i];
242 : }
243 : }
244 : };
245 :
246 : ///\ref levmarInterface fitter structure for the shifted Gamma Distribution
247 : /**
248 : *
249 : * \ingroup gammaDist_peak_fit
250 : *
251 : */
252 : template <typename _realT>
253 : struct gammaDistribution_3param_fitter
254 : {
255 : typedef _realT realT;
256 :
257 : static const int nparams = 3;
258 :
259 : static void func( realT *p, realT *hx, int m, int n, void *adata )
260 : {
261 : array2FitGammaDistribution<realT> *arr = (array2FitGammaDistribution<realT> *)adata;
262 :
263 : realT denom = func::gammaDistributionDenom<realT>( p[0], p[1] );
264 :
265 : for( int i = 0; i < arr->n; i++ )
266 : {
267 : hx[i] = func::gammaDistribution<realT>( i, p[2], p[0], p[1], denom ) - arr->data[i];
268 : }
269 : }
270 : };
271 :
272 : ///\ref levmarInterface fitter structure for the Gamma Distribution
273 : /**
274 : *
275 : * \ingroup gammaDist_peak_fit
276 : *
277 : */
278 : template <typename _realT>
279 : struct gammaDistribution_2param_fitter
280 : {
281 : typedef _realT realT;
282 :
283 : static const int nparams = 2;
284 :
285 : static void func( realT *p, realT *hx, int m, int n, void *adata )
286 : {
287 : array2FitGammaDistribution<realT> *arr = (array2FitGammaDistribution<realT> *)adata;
288 :
289 : realT denom = func::gammaDistributionDenom<realT>( p[0], p[1] );
290 :
291 : for( int i = 0; i < arr->n; i++ )
292 : {
293 : hx[i] = func::gammaDistribution<realT>( i, 0.0, p[0], p[1], denom ) - arr->data[i];
294 : }
295 : }
296 : };
297 :
298 : } // namespace fit
299 : } // namespace math
300 : } // namespace mx
301 :
302 : #endif // fitGammaDistribution_hpp
|
