normal.cpp

Go to the documentation of this file.

/**
 * @file normal.cpp
 * @brief Normal(μ, σ) family implementation.
 *
 * Self-contained family implementation: the class is file-local and
 * reaches the evaluators only through the registrars at the bottom
 * (DistributionRegistry descriptor + pairwise comparator / closure
 * rules).
 */
#include "DistributionCommon.h"
 
#include <algorithm>
#include <cmath>
#include <memory>
#include <optional>
#include <random>
#include <string>
#include <utility>
#include <vector>
 
namespace provsql {
 
namespace {
 
// (j-1)!! with the empty-product convention (-1)!! = 1.
double double_factorial_minus_one(unsigned j)
{
  if (j == 0) return 1.0;
  double r = 1.0;
  for (unsigned i = 1; i < j; i += 2) r *= static_cast<double>(i);
  return r;
}
 
/**
 * @brief Raw moments of @c X ~ Normal(μ, σ) truncated to @c [a, b].
 *
 * Closed form via the integration-by-parts recurrence on the
 * standardised variable Z = (X - μ)/σ:
 *   E[Z^{k}|α<Z<β] = (k-1) E[Z^{k-2}|α<Z<β]
 *                  + (α^{k-1}φ(α) − β^{k-1}φ(β)) / (Φ(β) − Φ(α))
 * with E[Z^0|…] = 1 and E[Z^1|…] = (φ(α) − φ(β)) / (Φ(β) − Φ(α))
 * (Greene, "Econometric Analysis", 5e, App. F).  Then expand
 * E[X^k] = E[(μ + σZ)^k] binomially.
 *
 * @c α = -∞ corresponds to @p a = -INFINITY (semi-infinite left tail);
 * @c β = +∞ to @p b = +INFINITY.  Returns @c NaN if @c P(α<Z<β) is
 * below a numerical floor (so the caller falls through to MC).
 */
double truncated_normal_raw_moment(double mu, double sigma, double a, double b,
                                   unsigned k)
{
  const double alpha = std::isfinite(a) ? (a - mu) / sigma : -kInf;
  const double beta  = std::isfinite(b) ? (b - mu) / sigma : +kInf;
  const double Phi_alpha = std::isfinite(alpha) ? Phi(alpha) : 0.0;
  const double Phi_beta  = std::isfinite(beta)  ? Phi(beta)  : 1.0;
  const double Z = Phi_beta - Phi_alpha;
  if (Z < 1e-12) return kNaN;
 
  const double phi_alpha = std::isfinite(alpha) ? phi(alpha) : 0.0;
  const double phi_beta  = std::isfinite(beta)  ? phi(beta)  : 0.0;
 
  /* E[Z^k | α<Z<β] via recurrence; store all moments up to k. */
  std::vector<double> M(k + 1, 0.0);
  M[0] = 1.0;
  if (k >= 1) M[1] = (phi_alpha - phi_beta) / Z;
  for (unsigned m = 2; m <= k; ++m) {
    /* α^{m-1}·φ(α) and β^{m-1}·φ(β); take 0 when the endpoint is
     * infinite (the φ factor vanishes faster than any polynomial). */
    double end_term = 0.0;
    if (std::isfinite(alpha))
      end_term += std::pow(alpha, static_cast<double>(m - 1)) * phi_alpha;
    if (std::isfinite(beta))
      end_term -= std::pow(beta, static_cast<double>(m - 1)) * phi_beta;
    M[m] = (m - 1) * M[m - 2] + end_term / Z;
  }
 
  /* E[X^k] = E[(μ + σZ)^k] = Σ_{i=0..k} C(k,i) μ^{k-i} σ^i E[Z^i|…]. */
  double total = 0.0;
  for (unsigned i = 0; i <= k; ++i) {
    total += binomial_coeff(k, i)
           * std::pow(mu, static_cast<double>(k - i))
           * std::pow(sigma, static_cast<double>(i))
           * M[i];
  }
  return total;
}
 
/** @brief Normal(μ=p1, σ=p2). */
class NormalDistribution final : public BaseDistribution {
public:
  using BaseDistribution::BaseDistribution;
  const DistributionFamily &family() const override;
  double mean() const override { return p1_; }
  bool meanIsAffine() const override { return true; }   // mean = μ
  double variance() const override { return p2_ * p2_; }
  double rawMoment(unsigned k) const override {
    if (k == 0) return 1.0;
    if (k == 1) return mean();
    // E[X^k] = sum_{j=0,2,...}^{k} C(k,j) mu^(k-j) sigma^j (j-1)!!
    const double mu = p1_, sigma = p2_;
    double total = 0.0;
    for (unsigned j = 0; j <= k; j += 2) {
      total += binomial_coeff(k, j)
             * std::pow(mu, static_cast<double>(k - j))
             * std::pow(sigma, static_cast<double>(j))
             * double_factorial_minus_one(j);
    }
    return total;
  }
  double pdf(double c) const override {
    if (!(p2_ > 0.0)) return kNaN;
    const double SQRT_2PI = std::sqrt(2.0 * M_PI);
    const double z = (c - p1_) / p2_;
    return std::exp(-0.5 * z * z) / (p2_ * SQRT_2PI);
  }
  double cdf(double c) const override {
    const double SQRT2 = std::sqrt(2.0);
    const double z = (c - p1_) / p2_;
    return 0.5 * (1.0 + std::erf(z / SQRT2));
  }
  DistSupport support() const override { return {-kInf, kInf}; }
  bool integrationRange(double &lo, double &hi) const override {
    if (!(p2_ > 0.0)) return false;
    lo = p1_ - 12.0 * p2_;
    hi = p1_ + 12.0 * p2_;
    return true;
  }
  std::pair<double, double> plotRange(double trunc_lo, double trunc_hi) const override {
    double lo = trunc_lo, hi = trunc_hi;
    if (!std::isfinite(lo)) lo = p1_ - 4.0 * p2_;
    if (!std::isfinite(hi)) hi = p1_ + 4.0 * p2_;
    return {lo, hi};
  }
  double sample(std::mt19937_64 &rng) const override {
    std::normal_distribution<double> d(p1_, p2_);
    return d(rng);
  }
  std::optional<double> quantile(double p) const override {
    if (!(p2_ > 0.0)) return std::nullopt;
    /* Beasley-Springer-Moro start (~1e-7), polished to machine
     * precision by two Newton steps on Φ(z) - p (φ(z) > 0 everywhere,
     * so the iteration is well-defined; the first step already lands
     * within ~1e-14). */
    double z = inv_phi(p);
    for (int i = 0; i < 2; ++i) {
      const double d = phi(z);
      if (!(d > 0.0)) break;   /* far tail: φ underflowed, keep BSM */
      z -= (Phi(z) - p) / d;
    }
    return p1_ + p2_ * z;
  }
  std::optional<double> truncatedRawMoment(double lo, double hi,
                                           unsigned k) const override {
    const double r = truncated_normal_raw_moment(p1_, p2_, lo, hi, k);
    if (std::isnan(r)) return std::nullopt;
    return r;
  }
  std::optional<std::vector<double>> sampleTruncated(
    std::mt19937_64 &rng, double lo, double hi, unsigned n) const override {
    const double mu = p1_, sigma = p2_;
    if (!(sigma > 0.0)) return std::nullopt;
    const double sqrt2 = std::sqrt(2.0);
    const double alpha = (lo - mu) / sigma;
    const double beta  = (hi - mu) / sigma;
    const double Phi_a = std::isfinite(alpha)
      ? 0.5 * (1.0 + std::erf(alpha / sqrt2))
      : (alpha < 0 ? 0.0 : 1.0);
    const double Phi_b = std::isfinite(beta)
      ? 0.5 * (1.0 + std::erf(beta / sqrt2))
      : (beta < 0 ? 0.0 : 1.0);
    if (!(Phi_a < Phi_b)) return std::nullopt;
    /* Inverse-CDF transform.  Clamp the target probability strictly
     * inside (0, 1) so @c inv_phi does not diverge near the asymptotes.
     * The 1e-15 margin is well below the BSM approximation's intrinsic
     * accuracy floor (~1e-7), so it's a safe sentinel. */
    static constexpr double EPS = 1e-15;
    std::uniform_real_distribution<double> U01(0.0, 1.0);
    std::vector<double> out;
    out.reserve(n);
    for (unsigned i = 0; i < n; ++i) {
      double u = Phi_a + U01(rng) * (Phi_b - Phi_a);
      if (u < EPS)        u = EPS;
      if (u > 1.0 - EPS)  u = 1.0 - EPS;
      const double z = inv_phi(u);
      out.push_back(mu + sigma * z);
    }
    return out;
  }
  std::unique_ptr<Distribution> affine(double a, double b) const override {
    if (a == 0.0) return nullptr;
    double mu = a * p1_;
    if (b != 0.0) mu += b;
    return std::make_unique<NormalDistribution>(mu, std::fabs(a) * p2_);
  }
  std::string serialise() const override {
    return "normal:" + double_to_text(p1_) + "," + double_to_text(p2_);
  }
  std::optional<double> asDirac() const override {
    if (p2_ == 0.0) return p1_;
    return std::nullopt;
  }
};
 
const DistributionFamily normal_family = {
  "normal", 2, "N", {"μ", "σ"},
  +[](double p1, double p2) -> std::unique_ptr<Distribution> {
    return std::make_unique<NormalDistribution>(p1, p2);
  }};
 
const DistributionFamily &NormalDistribution::family() const
{
  return normal_family;
}
 
/* (Normal, +, Normal): every term a·Zᵢ + bᵢ over independent normals
 * sums to Normal(Σ(aᵢμᵢ + bᵢ), √(Σ aᵢ²σᵢ²)).  A degenerate total
 * variance would make the closure a Dirac at the total mean; decline and
 * leave that to other passes (in practice unreachable: σ=0 normals are
 * constructed as gate_value by provsql.normal, so total_var > 0 whenever
 * the closure fires). */
std::unique_ptr<Distribution>
normalSumRule(const std::vector<ClosureTerm> &terms)
{
  double total_mean = 0.0;
  double total_var  = 0.0;
  for (const auto &t : terms) {
    total_mean += t.b;
    if (!t.dist) continue;
    const double mu    = t.dist->p1();
    const double sigma = t.dist->p2();
    total_mean += t.a * mu;
    total_var  += t.a * t.a * sigma * sigma;
  }
  if (total_var <= 0.0) return nullptr;
  return std::make_unique<NormalDistribution>(total_mean,
                                              std::sqrt(total_var));
}
 
[[maybe_unused]] const ClosureRuleRegistrar normal_sum_rule(
  "normal", "normal", &normalSumRule);
 
/* P(X < Y) for X, Y independent normals.  Reduces to P(X - Y < 0) with
 * X - Y ~ N(μ_X - μ_Y, √(σ_X² + σ_Y²)). */
double normalPairLess(const Distribution &X, const Distribution &Y)
{
  return NormalDistribution(
    X.p1() - Y.p1(),
    std::sqrt(X.p2() * X.p2() + Y.p2() * Y.p2())).cdf(0.0);
}
 
[[maybe_unused]] const ComparatorRuleRegistrar normal_less_rule(
  "normal", "normal", &normalPairLess);
 
/* Normal-Normal conjugacy: a Normal(θ, σ) observation with the latent in
 * the mean slot updates a Normal(μ, σ₀) prior by precision weighting:
 * τ ← 1/σ₀² + 1/σ², μ ← (μ/σ₀² + d/σ²)/τ.  The predictive of the datum
 * before the update is Normal(μ, √(σ₀² + σ²)) -- the marginal of
 * d = θ + noise. */
bool normalMeanConjugateUpdate(double &mu, double &sigma,
                               const DistributionTemplate &lik, double d)
{
  const double s = lik.p2.literal;
  if (!(s > 0.0) || !(sigma > 0.0)) return false;
  const double tau0 = 1.0 / (sigma * sigma);
  const double taul = 1.0 / (s * s);
  mu = (tau0 * mu + taul * d) / (tau0 + taul);
  sigma = std::sqrt(1.0 / (tau0 + taul));
  return true;
}
 
double normalMeanLogPredictive(double mu, double sigma,
                               const DistributionTemplate &lik, double d)
{
  const double s = lik.p2.literal;
  if (!(s > 0.0) || !(sigma > 0.0)) return kNaN;
  const double v = sigma * sigma + s * s;
  return -0.5 * ((d - mu) * (d - mu) / v + std::log(2.0 * M_PI * v));
}
 
[[maybe_unused]] const ConjugateRuleRegistrar normal_mean_conjugate(
  "normal", 0, "normal",
  {&normalMeanConjugateUpdate, &normalMeanLogPredictive});
 
[[maybe_unused]] const DistributionFamilyRegistrar normal_family_registrar(
  normal_family);
 
}  // namespace
 
}  // namespace provsql