bivariate_normal.hpp

Go to the documentation of this file.

#pragma once
 
#include <randlib/math/constants.hpp>
#include <randlib/math/functions.hpp>
#include <randlib/core/errors.hpp>
#include <randlib/core/continuous_bivariate.hpp>
#include <randlib/continuous/normal.hpp>
 
namespace randlib
{
 
template <typename RealType>
class BivariateNormalRand final
    : public ContinuousBivariateDistribution<NormalRand<RealType>, NormalRand<RealType>, RealType>
{
    Pair<RealType> variate(detail::EngineRef g) const override { return (*this)(g); }
 
    double mu1 = 0; 
    double mu2 = 0; 
    double sigma1 = 1; 
    double sigma2 = 1; 
    double rho = 0; 
    double pdfCoef = 1 + randlib::consts::ln2 + randlib::consts::lnpi; 
    double sqrt1mroSq = 1; 
 
public:
    BivariateNormalRand(double location1, double location2, double scale1, double scale2, double correlation);
    [[nodiscard]] std::string Name() const override;
 
    void SetLocations(double location1, double location2);
    void SetCovariance(double scale1, double scale2, double correlation);
    [[nodiscard]] inline DoublePair GetLocation() { return this->Mean(); }
    [[nodiscard]] inline double GetFirstLocation() const noexcept { return mu1; }
    [[nodiscard]] inline double GetSecondLocation() const noexcept { return mu2; }
    [[nodiscard]] inline double GetFirstScale() const noexcept { return sigma1; }
    [[nodiscard]] inline double GetSecondScale() const noexcept { return sigma2; }
    [[nodiscard]] inline double GetCorrelation() const noexcept { return rho; }
 
    [[nodiscard]] double f(const Pair<RealType> &point) const override;
    [[nodiscard]] double logf(const Pair<RealType> &point) const override;
    [[nodiscard]] double F(const Pair<RealType> &point) const override;
    template <std::uniform_random_bit_generator G>
    [[nodiscard]] Pair<RealType> operator()(G &g) const;
 
    long double Correlation() const override;
 
    [[nodiscard]] Pair<RealType> Mode() const override;
};
 
} // namespace randlib
 
namespace randlib
{
 
template <typename RealType>
BivariateNormalRand<RealType>::BivariateNormalRand(double location1, double location2, double scale1, double scale2,
                                                   double correlation)
{
    SetLocations(location1, location2);
    SetCovariance(scale1, scale2, correlation);
}
 
template <typename RealType>
std::string BivariateNormalRand<RealType>::Name() const
{
    return "Bivariate Normal( (" + this->toStringWithPrecision(GetFirstLocation()) + ", " +
           this->toStringWithPrecision(GetSecondLocation()) + "), (" + this->toStringWithPrecision(GetFirstScale()) +
           ", " + this->toStringWithPrecision(GetSecondScale()) + ", " + this->toStringWithPrecision(GetCorrelation()) +
           ") )";
}
 
template <typename RealType>
void BivariateNormalRand<RealType>::SetLocations(double location1, double location2)
{
    mu1 = location1;
    mu2 = location2;
    this->X.SetLocation(mu1);
    this->Y.SetLocation(mu2);
}
 
template <typename RealType>
void BivariateNormalRand<RealType>::SetCovariance(double scale1, double scale2, double correlation)
{
    if (!(scale1 > 0.0 && scale2 > 0.0))
        throw std::invalid_argument(
            "Scales of Bivariate-Normal distribution should be positive, but they're equal to " +
            std::to_string(scale1) + " and " + std::to_string(scale2) + " respectively");
    if (!(std::fabs(correlation) <= 1.0))
        throw std::invalid_argument(
            "Correlation of Bivariate-Normal distribution should be in interval [-1, 1], but it's equal to " +
            std::to_string(correlation));
    sigma1 = scale1;
    sigma2 = scale2;
    rho = correlation;
 
    this->X.SetScale(sigma1);
    this->Y.SetScale(sigma2);
 
    sqrt1mroSq = 0.5 * std::log1pl(-rho * rho);
    pdfCoef = this->X.GetLogScale() + this->Y.GetLogScale() + sqrt1mroSq + randlib::consts::ln2 + randlib::consts::lnpi;
    sqrt1mroSq = std::exp(sqrt1mroSq);
}
 
template <typename RealType>
double BivariateNormalRand<RealType>::f(const Pair<RealType> &point) const
{
    RealType xAdj = (point.first - mu1) / sigma1, yAdj = (point.second - mu2) / sigma2;
    if (rho == 1.0) 
        return (xAdj - yAdj == 0) ? INFINITY : 0.0;
    if (rho == -1.0) 
        return (xAdj + yAdj == 0) ? INFINITY : 0.0;
    return std::exp(logf(point));
}
 
template <typename RealType>
double BivariateNormalRand<RealType>::logf(const Pair<RealType> &point) const
{
    if (rho == 0.0) 
        return this->X.logf(point.first) + this->Y.logf(point.second);
    RealType xAdj = (point.first - mu1) / sigma1, yAdj = (point.second - mu2) / sigma2;
    if (rho == 1.0) 
        return (xAdj - yAdj == 0) ? INFINITY : -INFINITY;
    if (rho == -1.0) 
        return (xAdj + yAdj == 0) ? INFINITY : -INFINITY;
    RealType z = xAdj * xAdj - 2 * rho * xAdj * yAdj + yAdj * yAdj;
    return -(pdfCoef + 0.5 * z / (sqrt1mroSq * sqrt1mroSq));
}
 
template <typename RealType>
double BivariateNormalRand<RealType>::F(const Pair<RealType> &point) const
{
    if (rho == 0.0) 
        return this->X.F(point.first) * this->Y.F(point.second);
    double p1, p2 = 0.0;
    RealType x = point.first, y = point.second;
    RealType xAdj = (x - mu1) / sigma1, yAdj = (y - mu2) / sigma2;
    if (rho > 0)
        p1 = (xAdj < yAdj) ? this->X.F(x) : this->Y.F(y);
    else {
        double standardComplement = 0.5 * std::erfc(yAdj * randlib::consts::inv_sqrt2);
        p1 = std::max(this->X.F(x) - standardComplement, 0.0);
    }
    if (std::fabs(rho) < 1.0) {
        RealType lowLimit = std::asin(rho);
        RealType highLimit = math::sign(rho) * randlib::consts::half_pi;
        p2 = math::integrate(
            [xAdj, yAdj](double theta) {
                RealType cosTheta = std::cos(theta);
                RealType tanTheta = std::tan(theta);
                RealType integrand = xAdj * tanTheta - yAdj / cosTheta;
                integrand *= integrand;
                integrand += xAdj * xAdj;
                integrand = std::exp(-0.5 * integrand);
                return integrand;
            },
            lowLimit, highLimit);
    }
    return p1 - 0.5 * p2 / randlib::consts::pi;
}
 
template <typename RealType>
template <std::uniform_random_bit_generator G>
Pair<RealType> BivariateNormalRand<RealType>::operator()(G &g) const
{
    RealType Z1 = detail::standardNormal<RealType>(g);
    RealType Z2 = detail::standardNormal<RealType>(g);
    RealType x = mu1 + sigma1 * Z1;
    RealType y = mu2 + sigma2 * (rho * Z1 + sqrt1mroSq * Z2);
    return std::make_pair(x, y);
}
 
template <typename RealType>
long double BivariateNormalRand<RealType>::Correlation() const
{
    return rho;
}
 
template <typename RealType>
Pair<RealType> BivariateNormalRand<RealType>::Mode() const
{
    return std::make_pair(mu1, mu2);
}
 
} // namespace randlib