expModGaussian.hpp

Go to the documentation of this file.

/** \file expModGaussian.hpp
 * \brief The Exponentially Modified Gaussian distribution.
 * \ingroup gen_math_files
 * \author Jared R. Males (jaredmales@gmail.com)
 *
 */
 
//***********************************************************************//
// Copyright 2023 Jared R. Males (jaredmales@gmail.com)
//
// This file is part of mxlib.
//
// mxlib is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// mxlib is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with mxlib.  If not, see <http://www.gnu.org/licenses/>.
//***********************************************************************//
 
#ifndef expModGaussian_hpp
#define expModGaussian_hpp
 
#include <boost/math/tools/minima.hpp>
#include <limits>
 
namespace mx
{
namespace math
{
namespace func
{
 
/// The Exponentially Modified Gaussian at a point.
/** Calculates the value of the Exponentially Modified Gaussian distribution at a location specified by x.
 *
 *
 * \tparam realT a real floating point type
 *
 * \returns the value of the Exponentially Modified Gaussian distribution at x.
 *
 * \ingroup gen_math_expModGaussian
 */
template <typename realT>
realT expModGaussian( realT x,     ///< [in] the location at which to calculate the distribution
                      realT mu,    ///< [in] the mean parameter
                      realT sigma, ///< [in] the standard deviation
                      realT lambda ///< [in] the rate of decay
)
{
    return ( lambda / 2 ) * exp( ( lambda / 2 ) * ( 2 * mu + lambda * sigma * sigma - 2 * x ) ) *
           std::erfc( ( mu + lambda * sigma * sigma - x ) / ( root_two<realT>() * sigma ) );
}
 
/// The Mean of the Exponentially Modified Gaussian.
/** Calculates the mean of the Exponentially Modified Gaussian distribution.
 *
 *
 * \tparam realT a real floating point type
 *
 * \returns the mean of the Exponentially Modified Gaussian.
 *
 * \ingroup gen_math_expModGaussian
 */
template <typename realT>
realT expModGaussianMean( realT mu,    ///< [in] the mean parameter
                          realT lambda ///< [in] the rate of decay
)
{
    return mu + 1.0 / lambda;
}
 
/// The Variance of the Exponentially Modified Gaussian.
/** Calculates the variance of the Exponentially Modified Gaussian distribution.
 *
 *
 * \tparam realT a real floating point type
 *
 * \returns the variance of the Exponentially Modified Gaussian.
 *
 * \ingroup gen_math_expModGaussian
 */
template <typename realT>
realT expModGaussianVariance( realT sigma, ///< [in] the standard deviation
                              realT lambda ///< [in] the rate of decay
)
{
    return sigma * sigma + 1.0 / ( lambda * lambda );
}
 
template <typename realT>
struct emgModeFunc
{
    realT mu;
    realT sigma;
    realT lambda;
 
    realT operator()( const realT &x )
    {
        return -expModGaussian( x, mu, sigma, lambda );
    }
};
 
/// The Mode of the Exponentially Modified Gaussian.
/** Calculates the mode of the Exponentially Modified Gaussian distribution.
 * This is done iteratively with Brent's method.
 *
 * \tparam realT a real floating point type
 *
 * \returns the mode of the Exponentially Modified Gaussian.
 *
 * \ingroup gen_math_expModGaussian
 */
template <typename realT>
realT expModGaussianMode( realT mu,    ///< [in] the mean parameter
                          realT sigma, ///< [in] the standard deviation
                          realT lambda ///< [in] the rate of decay
)
{
    realT mn = expModGaussianMean( mu, lambda );
    realT sd = sqrt( expModGaussianVariance( sigma, lambda ) );
 
    std::cerr << mn << " " << sd << "\n";
 
    emgModeFunc<realT> mf;
    mf.mu = mu;
    mf.sigma = sigma;
    mf.lambda = lambda;
 
    uintmax_t maxit = 1000;
    try
    {
        std::pair<realT, realT> brack;
        brack = boost::math::tools::brent_find_minima<emgModeFunc<realT>, realT>(
            mf, mn - 2 * sd, mn + 2 * sd, std::numeric_limits<realT>::digits, maxit );
        std::cerr << brack.first << " " << brack.second << " " << maxit << "\n";
 
        return brack.first;
    }
    catch( ... )
    {
        std::cerr << "expModGaussianMode: No mode found\n";
        return std::numeric_limits<realT>::quiet_NaN();
    }
}
 
} // namespace func
} // namespace math
} // namespace mx
 
#endif // expModGaussian_hpp