normal distribution

ソースコード

 1#include <iostream>
 2#include <cmath>
 3#include <cassert>
 4using namespace std;
 5
 6namespace titan23 {
 7
 8class NormalDist {
 9private:
10    static constexpr const double sqrt2 = 1.4142135623730951;
11    static constexpr const double a1 = 0.254829592, a2 = -0.284496736, a3 = 1.421413741, a4 = -1.453152027, a5 = 1.061405429;
12    static constexpr const double p1 = 0.3275911;
13
14    static double cdf_approx(double x) {
15        double sgn = (x < 0) ? -1.0 : 1.0;
16        x = std::fabs(x) / std::sqrt(2.0);
17        double t = 1.0 / (1.0 + p1 * x);
18        double y = 1.0 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * std::exp(-x * x);
19        return 0.5 * (1.0 + sgn * y);
20    }
21
22    static double erfinv(double x) {
23        /// 逆誤差関数 erfinv の近似計算
24        assert(-1.0 <= x && x <= 1.0);
25        double a = 0.0;
26        double b = 0.0;
27        double t = x;
28        double y = t;
29        double epsilon = 1e-10;
30        while (true) {
31            t = (2.0 * y + x / (sqrt2 * std::exp(-y * y))) / 2.0;
32            if (std::fabs(t - y) < epsilon) break;
33            y = t;
34        }
35        return y;
36    }
37
38public:
39    /// 正規分布 N(mu, sigma) において、x 以下となる確率を返す
40    static double cdf(double x, double mu, double sigma) {
41        return 0.5 * (1 + std::erf((x - mu) / (sigma * sqrt2)));
42    }
43
44    /// 正規分布 N(mu, sigma) において、x 以下となる確率を返す
45    static double cdf_approx(double x, double mu, double sigma) {
46        return cdf_approx((x - mu) / sigma);
47    }
48
49    /// 正規分布 N(mu, sigma) における x の確率密度を返す
50    static double pdf(double x, double mu, double sigma) {
51        double coeff = 1.0 / (sigma * std::sqrt(2.0 * M_PI));
52        double exponent = -0.5 * std::pow((x - mu) / sigma, 2);
53        return coeff * std::exp(exponent);
54    }
55
56    /// 正規分布 N(mu, sigma) に対して、観測した値が区間 [a, b] に含まれる確率を返す
57    static double probability_in_range(double l, double r, double mu, double sigma) {
58        assert(l <= r);
59        assert(cdf(r, mu, sigma) - cdf(l, mu, sigma) >= 0);
60        return cdf(r, mu, sigma) - cdf(l, mu, sigma);
61    }
62
63    /// 正規分布 N(mu, sigma) における累積確率 p に対応する x を返す
64    static double inverse_cdf(double p, double mu, double sigma) {
65        assert(0 <= p && p <= 1.0);
66        return mu + sigma * sqrt2 * erfinv(2.0 * p - 1.0);
67    }
68};
69} // namespace titan23

仕様

Warning

doxygenfile: Cannot find file “titan_cpplib/ahc/normal_distribution.cpp