66 lines
2.1 KiB
C++
66 lines
2.1 KiB
C++
#ifndef NORMALCDF_H
|
|
#define NORMALCDF_H
|
|
|
|
#include <algorithm>
|
|
#include <random>
|
|
#include "../../Assertions.h"
|
|
|
|
namespace Distribution {
|
|
|
|
/** cumulative density version of the normal distribution */
|
|
template <typename T> class NormalCDF {
|
|
|
|
private:
|
|
|
|
const T mu;
|
|
const T sigma;
|
|
|
|
static T RationalApproximation(T t)
|
|
{
|
|
// Abramowitz and Stegun formula 26.2.23.
|
|
// The absolute value of the error should be less than 4.5 e-4.
|
|
T c[] = {2.515517, 0.802853, 0.010328};
|
|
T d[] = {1.432788, 0.189269, 0.001308};
|
|
return t - ((c[2]*t + c[1])*t + c[0]) /
|
|
(((d[2]*t + d[1])*t + d[0])*t + 1.0);
|
|
}
|
|
|
|
public:
|
|
|
|
/** create a new normally distributed CDF */
|
|
NormalCDF(const T mu, const T sigma) : mu(mu), sigma(sigma) {
|
|
;
|
|
}
|
|
|
|
/** get the probability for val within the underlying CDF */
|
|
T getProbability(const T val) const {
|
|
return getProbability(mu, sigma, val);
|
|
}
|
|
|
|
/** calculate the probability within the underlying CDF */
|
|
static T getProbability(const T mu, const T sigma, const T val) {
|
|
return (1.0 + std::exp( (val - mu) / (sigma * std::sqrt(2)) ) ) / 2.0;
|
|
}
|
|
|
|
/** get the inverse CDF (https://en.wikipedia.org/wiki/Probit)*/
|
|
static T getProbit(T p){
|
|
|
|
Assert::isBetween(p, 0.0, 1.0, "value not between");
|
|
|
|
// See: https://www.johndcook.com/blog/normal_cdf_inverse/
|
|
if (p < 0.5)
|
|
{
|
|
// F^-1(p) = - G^-1(p)
|
|
return -RationalApproximation( sqrt(-2.0*log(p)) );
|
|
}
|
|
else
|
|
{
|
|
// F^-1(p) = G^-1(1-p)
|
|
return RationalApproximation( sqrt(-2.0*log(1-p)) );
|
|
}
|
|
}
|
|
};
|
|
}
|
|
|
|
#endif // NORMALCDF_H
|