f-Divergences

1
#include <iostream>
2
#include <vector>
3
#include <cmath>
4
#include <functional>
5
#include <limits>
6
#include <numeric>
7
#include <stdexcept>
8
9
// Utility: check sums to ~1
10
bool approx_one(const std::vector<double>& v, double tol = 1e-9) {
11
    double s = std::accumulate(v.begin(), v.end(), 0.0);
12
    return std::fabs(s - 1.0) <= tol;
13
}
14
15
// Common f-generators as functions on [0, +inf)
16
namespace fgen {
17
    // KL: f(t) = t * log(t), with f(1)=0, limit f(0)=0
18
    double kl(double t) {
19
        if (t <= 0.0) return 0.0; // by limit t log t -> 0 as t -> 0+
20
        return t * std::log(t);
21
    }
22
    // Total Variation: f(t) = 0.5 * |t - 1|
23
    double tv(double t) {
24
        return 0.5 * std::fabs(t - 1.0);
25
    }
26
    // Hellinger squared: f(t) = (sqrt(t) - 1)^2, with f(0)=1
27
    double hellinger2(double t) {
28
        if (t <= 0.0) return 1.0; // limit as t->0+: (0 - 1)^2 = 1
29
        double s = std::sqrt(t);
30
        return (s - 1.0) * (s - 1.0);
31
    }
32
    // Pearson chi-square: f(t) = (t - 1)^2
33
    double chi2(double t) {
34
        double d = t - 1.0;
35
        return d * d;
36
    }
37
    // Jensen-Shannon (base e): f(t) = t log t - (t+1) log((t+1)/2) + log 2
38
    double js(double t) {
39
        if (t <= 0.0) return std::log(2.0); // limit t->0+: 0 - (1)log(1/2) + log 2 = log 2
40
        double a = t * std::log(t);
41
        double b = (t + 1.0) * std::log((t + 1.0) / 2.0);
42
        return a - b + std::log(2.0);
43
    }
44
}
45
46
// Compute D_f(P||Q) for discrete distributions (finite support)
47
// P and Q must be the same size, entries >= 0, sum to 1 (approximately), and P << Q.
48
// If exists i with Q[i]==0 and P[i]>0, returns +inf.
49
double f_divergence_discrete(
50
    const std::vector<double>& P,
51
    const std::vector<double>& Q,
52
    const std::function<double(double)>& f
53
) {
54
    if (P.size() != Q.size() || P.empty()) {
55
        throw std::invalid_argument("P and Q must be non-empty and of the same size");
56
    }
57
    if (!approx_one(P) || !approx_one(Q)) {
58
        std::cerr << "Warning: P or Q does not sum to 1 within tolerance.\n";
59
    }
60
    double D = 0.0;
61
    for (size_t i = 0; i < P.size(); ++i) {
62
        double p = P[i], q = Q[i];
63
        if (q == 0.0) {
64
            if (p == 0.0) {
65
                continue; // contributes 0
66
            } else {
67
                return std::numeric_limits<double>::infinity(); // support mismatch
68
            }
69
        }
70
        double t = (p == 0.0) ? 0.0 : (p / q);
71
        D += q * f(t);
72
    }
73
    return D;
74
}
75
76
int main() {
77
    // Example: two categorical distributions over 4 outcomes
78
    std::vector<double> P = {0.10, 0.20, 0.30, 0.40};
79
    std::vector<double> Q = {0.25, 0.25, 0.25, 0.25}; // uniform
80
81
    std::cout.setf(std::ios::fixed); std::cout.precision(6);
82
83
    double D_KL = f_divergence_discrete(P, Q, fgen::kl);
84
    double D_TV = f_divergence_discrete(P, Q, fgen::tv);
85
    double D_H2 = f_divergence_discrete(P, Q, fgen::hellinger2);
86
    double D_CHI2 = f_divergence_discrete(P, Q, fgen::chi2);
87
    double D_JS = f_divergence_discrete(P, Q, fgen::js);
88
89
    std::cout << "KL(P||Q)           = " << D_KL  << "\n";
90
    std::cout << "TV(P,Q)            = " << D_TV  << "\n";
91
    std::cout << "Hellinger^2(P,Q)   = " << D_H2  << "\n";
92
    std::cout << "Pearson chi^2      = " << D_CHI2<< "\n";
93
    std::cout << "Jensen-Shannon     = " << D_JS  << "\n";
94
95
    return 0;
96
}
97

1
#include <iostream>
2
#include <random>
3
#include <cmath>
4
#include <functional>
5
#include <limits>
6
7
struct Normal {
8
    double mu;
9
    double sigma; // > 0
10
};
11
12
double log_pdf_normal(const Normal& dist, double x) {
13
    const double var = dist.sigma * dist.sigma;
14
    return -0.5 * std::log(2.0 * M_PI * var) - 0.5 * (x - dist.mu) * (x - dist.mu) / var;
15
}
16
17
// Safe exponential to avoid overflow/underflow extremes
18
inline double safe_exp(double a) {
19
    if (a > 700.0) return std::numeric_limits<double>::infinity();
20
    if (a < -745.0) return 0.0;
21
    return std::exp(a);
22
}
23
24
namespace fgen {
25
    double kl(double t) { return (t <= 0.0) ? 0.0 : t * std::log(t); }
26
    double hellinger2_from_logw(double logw) {
27
        // f(t) = (sqrt(t) - 1)^2 with t = exp(logw) => sqrt(t) = exp(0.5*logw)
28
        double s = safe_exp(0.5 * logw);
29
        double d = s - 1.0;
30
        return d * d;
31
    }
32
}
33
34
// Monte Carlo estimator using samples from Q: D_f(P||Q) ~ 1/n sum f(p(x)/q(x))
35
double estimate_f_divergence_mc(
36
    size_t n,
37
    const Normal& P,
38
    const Normal& Q,
39
    const std::function<double(double)>& f_from_t,
40
    std::mt19937& rng
41
) {
42
    std::normal_distribution<double> q_sampler(Q.mu, Q.sigma);
43
    double acc = 0.0;
44
    for (size_t i = 0; i < n; ++i) {
45
        double x = q_sampler(rng);
46
        double logw = log_pdf_normal(P, x) - log_pdf_normal(Q, x); // log(p/q)
47
        double t = safe_exp(logw);
48
        acc += f_from_t(t);
49
    }
50
    return acc / static_cast<double>(n);
51
}
52
53
// Closed-form KL between univariate Gaussians
54
double kl_gaussians_closed(const Normal& P, const Normal& Q) {
55
    double s1 = P.sigma, s2 = Q.sigma;
56
    double term = std::log(s2 / s1) + (s1*s1 + (P.mu - Q.mu)*(P.mu - Q.mu)) / (2.0 * s2 * s2) - 0.5;
57
    return term;
58
}
59
60
// Closed-form squared Hellinger between univariate Gaussians
61
double hellinger2_gaussians_closed(const Normal& P, const Normal& Q) {
62
    double s1 = P.sigma, s2 = Q.sigma;
63
    double num = 2.0 * s1 * s2;
64
    double den = s1*s1 + s2*s2;
65
    double expo = std::exp(- (P.mu - Q.mu)*(P.mu - Q.mu) / (4.0 * den));
66
    return 1.0 - std::sqrt(num / den) * expo;
67
}
68
69
int main() {
70
    std::mt19937 rng(123);
71
    Normal P{0.0, 1.0};
72
    Normal Q{1.0, 2.0};
73
74
    size_t n = 200000; // number of MC samples
75
76
    // Estimate KL via Monte Carlo using f(t) = t log t
77
    double Dkl_mc = estimate_f_divergence_mc(n, P, Q, [](double t){ return (t<=0.0)?0.0: t*std::log(t); }, rng);
78
    double Dkl_cf = kl_gaussians_closed(P, Q);
79
80
    // Estimate Hellinger^2 via Monte Carlo more stably using logw directly
81
    auto hell_from_t = [](double t){
82
        // fallback generic path if only t is available
83
        double s = std::sqrt(t);
84
        double d = s - 1.0; return d*d;
85
    };
86
    double Dh2_mc = estimate_f_divergence_mc(n, P, Q, hell_from_t, rng);
87
    double Dh2_cf = hellinger2_gaussians_closed(P, Q);
88
89
    std::cout.setf(std::ios::fixed); std::cout.precision(6);
90
    std::cout << "KL MC estimate      : " << Dkl_mc << "\n";
91
    std::cout << "KL closed-form      : " << Dkl_cf << "\n";
92
    std::cout << "Hellinger^2 MC est. : " << Dh2_mc << "\n";
93
    std::cout << "Hellinger^2 closed  : " << Dh2_cf << "\n";
94
95
    return 0;
96
}
97

1
#include <iostream>
2
#include <random>
3
#include <cmath>
4
#include <functional>
5
#include <limits>
6
7
struct Normal { double mu, sigma; };
8
9
double log_pdf_normal(const Normal& d, double x) {
10
    double v = d.sigma * d.sigma;
11
    return -0.5*std::log(2*M_PI*v) - 0.5*(x-d.mu)*(x-d.mu)/v;
12
}
13
14
inline double safe_exp(double a) {
15
    if (a > 700.0) return std::numeric_limits<double>::infinity();
16
    if (a < -745.0) return 0.0;
17
    return std::exp(a);
18
}
19
20
// Change-of-measure identity: E_Q[h] = E_P[h / w], where w = p/q
21
// Here h = f(w). So D_f(P||Q) = E_P[f(w)/w]
22
double estimate_via_change_of_measure(
23
    size_t n,
24
    const Normal& P,
25
    const Normal& Q,
26
    const std::function<double(double)>& f_from_t,
27
    std::mt19937& rng
28
) {
29
    std::normal_distribution<double> p_sampler(P.mu, P.sigma);
30
    double acc = 0.0;
31
    size_t finite = 0;
32
    for (size_t i = 0; i < n; ++i) {
33
        double x = p_sampler(rng);
34
        double logw = log_pdf_normal(P, x) - log_pdf_normal(Q, x); // log(p/q)
35
        double w = safe_exp(logw);
36
        if (w == 0.0 || !std::isfinite(w)) continue; // skip pathological contributions
37
        acc += f_from_t(w) / w; // unbiased if E[.|finite] well-defined
38
        ++finite;
39
    }
40
    return (finite == 0) ? std::numeric_limits<double>::quiet_NaN() : acc / static_cast<double>(finite);
41
}
42
43
int main(){
44
    std::mt19937 rng(42);
45
    Normal P{0.0, 1.0};
46
    Normal Q{1.0, 2.0};
47
48
    auto f_kl = [](double t){ return (t<=0.0)?0.0: t*std::log(t); };
49
50
    size_t n = 200000;
51
    double D_est = estimate_via_change_of_measure(n, P, Q, f_kl, rng);
52
53
    // Closed-form KL for comparison
54
    auto kl_cf = [](Normal P, Normal Q){
55
        return std::log(Q.sigma/P.sigma) + (P.sigma*P.sigma + (P.mu-Q.mu)*(P.mu-Q.mu))/(2*Q.sigma*Q.sigma) - 0.5;
56
    };
57
58
    std::cout.setf(std::ios::fixed); std::cout.precision(6);
59
    std::cout << "KL via change-of-measure (from P): " << D_est << "\n";
60
    std::cout << "KL closed-form                 : " << kl_cf(P,Q) << "\n";
61
62
    std::cout << "Note: This estimator can have high variance when w is small.\n";
63
    return 0;
64
}
65