Line data Source code
1 : /** \file expModGaussian.hpp
2 : * \brief The Exponentially Modified Gaussian distribution.
3 : * \ingroup gen_math_files
4 : * \author Jared R. Males (jaredmales@gmail.com)
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 expModGaussian_hpp
28 : #define expModGaussian_hpp
29 :
30 : #include <boost/math/tools/minima.hpp>
31 : #include <iostream>
32 : #include <limits>
33 :
34 : #include "../constants.hpp"
35 :
36 : namespace mx
37 : {
38 : namespace math
39 : {
40 : namespace func
41 : {
42 :
43 : /// The Exponentially Modified Gaussian at a point.
44 : /** Calculates the value of the Exponentially Modified Gaussian distribution at a location specified by x.
45 : *
46 : *
47 : * \tparam realT a real floating point type
48 : *
49 : * \returns the value of the Exponentially Modified Gaussian distribution at x.
50 : *
51 : * \ingroup gen_math_expModGaussian
52 : */
53 : template <typename realT>
54 0 : realT expModGaussian( realT x, ///< [in] the location at which to calculate the distribution
55 : realT mu, ///< [in] the mean parameter
56 : realT sigma, ///< [in] the standard deviation
57 : realT lambda ///< [in] the rate of decay
58 : )
59 : {
60 0 : return ( lambda / 2 ) * exp( ( lambda / 2 ) * ( 2 * mu + lambda * sigma * sigma - 2 * x ) ) *
61 0 : std::erfc( ( mu + lambda * sigma * sigma - x ) / ( root_two<realT>() * sigma ) );
62 : }
63 :
64 : /// The Mean of the Exponentially Modified Gaussian.
65 : /** Calculates the mean of the Exponentially Modified Gaussian distribution.
66 : *
67 : *
68 : * \tparam realT a real floating point type
69 : *
70 : * \returns the mean of the Exponentially Modified Gaussian.
71 : *
72 : * \ingroup gen_math_expModGaussian
73 : */
74 : template <typename realT>
75 : realT expModGaussianMean( realT mu, ///< [in] the mean parameter
76 : realT lambda ///< [in] the rate of decay
77 : )
78 : {
79 : return mu + 1.0 / lambda;
80 : }
81 :
82 : /// The Variance of the Exponentially Modified Gaussian.
83 : /** Calculates the variance of the Exponentially Modified Gaussian distribution.
84 : *
85 : *
86 : * \tparam realT a real floating point type
87 : *
88 : * \returns the variance of the Exponentially Modified Gaussian.
89 : *
90 : * \ingroup gen_math_expModGaussian
91 : */
92 : template <typename realT>
93 : realT expModGaussianVariance( realT sigma, ///< [in] the standard deviation
94 : realT lambda ///< [in] the rate of decay
95 : )
96 : {
97 : return sigma * sigma + 1.0 / ( lambda * lambda );
98 : }
99 :
100 : template <typename realT>
101 : struct emgModeFunc
102 : {
103 : realT mu;
104 : realT sigma;
105 : realT lambda;
106 :
107 0 : realT operator()( const realT &x )
108 : {
109 0 : return -expModGaussian( x, mu, sigma, lambda );
110 : }
111 : };
112 :
113 : /// The Mode of the Exponentially Modified Gaussian.
114 : /** Calculates the mode of the Exponentially Modified Gaussian distribution.
115 : * This is done iteratively with Brent's method.
116 : *
117 : * \tparam realT a real floating point type
118 : *
119 : * \returns the mode of the Exponentially Modified Gaussian.
120 : *
121 : * \ingroup gen_math_expModGaussian
122 : */
123 : template <typename realT>
124 : realT expModGaussianMode( realT mu, ///< [in] the mean parameter
125 : realT sigma, ///< [in] the standard deviation
126 : realT lambda ///< [in] the rate of decay
127 : )
128 : {
129 : realT mn = expModGaussianMean( mu, lambda );
130 : realT sd = sqrt( expModGaussianVariance( sigma, lambda ) );
131 :
132 : std::cerr << mn << " " << sd << "\n";
133 :
134 : emgModeFunc<realT> mf;
135 : mf.mu = mu;
136 : mf.sigma = sigma;
137 : mf.lambda = lambda;
138 :
139 : uintmax_t maxit = 1000;
140 : try
141 : {
142 : std::pair<realT, realT> brack;
143 : brack = boost::math::tools::brent_find_minima<emgModeFunc<realT>, realT>(
144 : mf, mn - 2 * sd, mn + 2 * sd, std::numeric_limits<realT>::digits, maxit );
145 : std::cerr << brack.first << " " << brack.second << " " << maxit << "\n";
146 :
147 : return brack.first;
148 : }
149 : catch( ... )
150 : {
151 : std::cerr << "expModGaussianMode: No mode found\n";
152 : return std::numeric_limits<realT>::quiet_NaN();
153 : }
154 : }
155 :
156 : } // namespace func
157 : } // namespace math
158 : } // namespace mx
159 :
160 : #endif // expModGaussian_hpp
|
