prior.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_LANG_PRIOR_H
#define LINE_LANG_PRIOR_H
 
/**
 * @file
 * @ingroup line_lang
 * `Prior`: parameter uncertainty as a weighted set of alternative models.
 *
 * Port of `matlab/src/lang/processes/Prior.m`. A Prior is placed where a
 * service or arrival distribution goes and says that the law is not known: it
 * is one of an explicit alternative set with prior weights (the DISCRETE form),
 * or it is generated by a scalar parameter whose density is given together with
 * a map from the parameter to a distribution (the CONTINUOUS form of Trivedi
 * and Bobbio (2017), Sec. 3.4). SolverUQ is the only consumer.
 *
 * IT IS NOT A MIXTURE, and the difference is not cosmetic. A mixture says each
 * JOB draws its service law afresh; a Prior says the MODEL has one law and we
 * do not know which. The two give different queue lengths, and it is the second
 * that epistemic uncertainty propagation means. The mixture moments are still
 * exposed (`prior_mean`, `prior_scv`, `prior_skewness`, `prior_cdf`,
 * `prior_lst`), because the reference's `Prior` inherits `Distribution` and
 * answers those questions, and because `refresh_rates` will lower a Prior to
 * SOME rate if a struct is dumped before SolverUQ has expanded it -- but no
 * solver ever reads them: `Feature::Prior` is declared by SolverUQ alone, so
 * every other solver refuses the model at the gate.
 *
 * DISCRETIZATION, which is where the two forms meet. `prior_discretize` reduces
 * either form to (dists, weights):
 *   quadrature   the discrete form UNCHANGED (it is already exact, and `n` is
 *                ignored); the continuous form at the conditional medians of n
 *                equal-mass strata, weights 1/n. The nodes are quantiles, so
 *                the rule integrates the parameter density in PROBABILITY space
 *                and needs only `dist_cdf`, which every family has.
 *   montecarlo   n draws, weights 1/n. For a discrete prior the draw is of the
 *                alternative INDEX against its probabilities -- returning the
 *                alternatives unweighted here would silently drop the prior.
 *
 * WHERE THE RANDOMNESS COMES FROM, and how it differs from the reference.
 * MATLAB draws with `rand` and, for the continuous form, with
 * `paramDist.sample(n)`, i.e. each family's own sampler. This port has no
 * per-family sampler at all -- `Distrib` carries none, and SSA generates its
 * clocks internally -- so the Monte Carlo design draws a uniform from an
 * explicit MT19937 stream and inverts the CDF. The design points are therefore
 * NOT the reference's for the same seed, exactly as the SSA sample paths are
 * not; the design's LAW is the same, and a quadrature design is identical
 * across the two codebases because it draws nothing.
 */
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <functional>
#include <memory>
#include <random>
#include <string>
#include <vector>
 
#include "line/lang/distribution.h"
#include "line/lang/lang_types.h"
#include "line/num/number.h"
#include "line/util/error.h"
 
namespace line {
namespace lang {
 
/** The default number of nodes per continuous Prior, MATLAB `options.samples`. */
inline constexpr std::size_t kPriorDefaultNodes = 11;
 
/**
 * The uniform stream the Monte Carlo design draws from.
 *
 * The same 53-bit assembly from two MT19937 words that `ssa::SsaRng` uses, so
 * the two random paths in this port share one construction; it is duplicated
 * rather than shared because `lang/` must not depend on a solver.
 */
class PriorRng {
  public:
    explicit PriorRng(unsigned long seed) : g_(static_cast<std::uint_fast32_t>(seed)) {}
    double uniform() {
        const std::uint64_t a = g_() >> 5, b = g_() >> 6;
        return (static_cast<double>(a) * 67108864.0 + static_cast<double>(b) + 0.5) /
               9007199254740992.0;
    }
 
  private:
    std::mt19937 g_;
};
 
/** A Prior reduced to alternatives and weights; the output of `prior_discretize`. */
template <class T>
struct PriorDesign {
    std::vector<Distrib<T>> dists;
    std::vector<T> weights;
};
 
/**
 * `Prior(distributions, probabilities)`: the discrete form.
 *
 * The weights must be nonnegative and sum to one within CoarseTol, which is the
 * reference's own tolerance; a set that does not is a modelling error and is
 * refused rather than renormalized, since renormalizing would answer for a
 * prior the caller did not write.
 */
template <class T>
Distrib<T> prior_discrete(const std::vector<Distrib<T>>& alternatives,
                          const std::vector<T>& probabilities) {
    const T zero = num_traits<T>::from_int(0), one = num_traits<T>::from_int(1);
    if (alternatives.empty()) throw InputError("Prior: the alternative set is empty");
    if (alternatives.size() != probabilities.size())
        throw InputError("Prior: there are " + std::to_string(alternatives.size()) +
                         " alternatives and " + std::to_string(probabilities.size()) +
                         " probabilities");
    T tot = zero;
    for (std::size_t i = 0; i < probabilities.size(); ++i) {
        if (probabilities[i] < zero) throw InputError("Prior: the probabilities must be nonnegative");
        tot += probabilities[i];
        if (alternatives[i].disabled)
            throw InputError("Prior: alternative " + std::to_string(i + 1) + " is disabled");
        if (alternatives[i].is_prior())
            throw InputError("Prior: alternative " + std::to_string(i + 1) +
                             " is itself a Prior; nest the uncertainty in one Prior instead");
    }
    if (std::fabs(num_traits<T>::to_double(T(tot - one))) > GlobalConstants::CoarseTol)
        throw InputError("Prior: the probabilities must sum to 1 (they sum to " +
                         std::to_string(num_traits<T>::to_double(tot)) + ")");
    Distrib<T> d;
    d.type = ProcessType::PRIOR;
    d.disabled = false;
    d.prior = std::make_shared<PriorSpec<T>>();
    d.prior->continuous = false;
    d.prior->alternatives = alternatives;
    d.prior->probabilities = probabilities;
    return d;
}
 
/**
 * `Prior(paramDist, distFactory)`: the continuous form.
 *
 * `factory` maps a value of the parameter to a distribution, e.g. a rate to an
 * `Exp`. It is called eagerly once here, on the median of the parameter law, so
 * that a factory returning a disabled or nested-Prior distribution is refused
 * at construction rather than at the first design point.
 */
template <class T>
Distrib<T> prior_continuous(const Distrib<T>& param_dist,
                            const std::function<Distrib<T>(const T&)>& factory) {
    if (param_dist.disabled) throw InputError("Prior: the parameter distribution is disabled");
    if (param_dist.is_prior())
        throw InputError("Prior: the parameter distribution is itself a Prior");
    if (!factory) throw InputError("Prior: the distribution factory is empty");
    const Distrib<T> probe = factory(dist_quantile(param_dist, num_traits<T>::from_double(0.5)));
    if (probe.disabled || probe.is_prior())
        throw InputError("Prior: the distribution factory must return a usable Distribution");
    Distrib<T> d;
    d.type = ProcessType::PRIOR;
    d.disabled = false;
    d.prior = std::make_shared<PriorSpec<T>>();
    d.prior->continuous = true;
    d.prior->param_dist = param_dist;
    d.prior->factory = factory;
    return d;
}
 
/**
 * `Prior.fromSample(k, s)`: the posterior of a rate estimated from lifetime data.
 *
 * Given k i.i.d. exponential observations summing to s, the Jeffreys prior
 * f(lambda) = s/lambda yields the posterior lambda^(k-1) s^k exp(-lambda s) /
 * (k-1)!, an Erlang density of k phases and phase rate s (Trivedi and Bobbio
 * (2017), Eq. 3.71). Its mean k/s is the maximum-likelihood rate and its
 * variance k/s^2 shrinks as k grows, so the prior concentrates on the estimate.
 */
template <class T>
Distrib<T> prior_from_sample(std::size_t k, const T& s,
                             const std::function<Distrib<T>(const T&)>& factory =
                                 std::function<Distrib<T>(const T&)>()) {
    if (k < 1) throw InputError("Prior.fromSample: k must be a positive observation count");
    if (!(s > num_traits<T>::from_int(0)))
        throw InputError("Prior.fromSample: s must be a positive sum of observed lifetimes");
    std::function<Distrib<T>(const T&)> f = factory;
    if (!f) f = [](const T& lambda) { return Distrib<T>::exp_rate(lambda); };
    return prior_continuous<T>(Distrib<T>::erlang(s, k), f);
}
 
/**
 * Reduce a Prior to `n` weighted alternatives, MATLAB `Prior.discretize`.
 *
 * `rng` is consulted only by the `montecarlo` method; a quadrature design draws
 * nothing and is reproducible across codebases.
 */
template <class T>
PriorDesign<T> prior_discretize(const Distrib<T>& d, std::size_t n, const std::string& method,
                                PriorRng& rng) {
    if (!d.is_prior()) throw InputError("prior_discretize: the distribution is not a Prior");
    if (method != "quadrature" && method != "montecarlo")
        throw InputError("prior_discretize: unknown discretization method '" + method + "'");
    if (n < 1) throw InputError("prior_discretize: the node count must be at least 1");
    const PriorSpec<T>& sp = *d.prior;
    PriorDesign<T> out;
    const T wn = T(num_traits<T>::from_int(1) / num_traits<T>::from_int(static_cast<long>(n)));
 
    if (!sp.continuous) {
        if (method == "quadrature") {
            out.dists = sp.alternatives;
            out.weights = sp.probabilities;
            return out;
        }
        // The draw is of the alternative INDEX against its own probabilities.
        std::vector<double> cum(sp.probabilities.size(), 0.0);
        double acc = 0.0;
        for (std::size_t i = 0; i < sp.probabilities.size(); ++i) {
            acc += num_traits<T>::to_double(sp.probabilities[i]);
            cum[i] = acc;
        }
        for (std::size_t i = 0; i < n; ++i) {
            const double u = rng.uniform() * acc;
            std::size_t k = 0;
            while (k + 1 < cum.size() && cum[k] < u) ++k;
            out.dists.push_back(sp.alternatives[k]);
            out.weights.push_back(wn);
        }
        return out;
    }
 
    for (std::size_t i = 0; i < n; ++i) {
        T p;
        if (method == "quadrature") {
            // The midpoint of stratum i, ((i+1) - 0.5)/n in the reference's
            // 1-based indexing: the conditional median of an equal-mass stratum.
            p = num_traits<T>::from_double((static_cast<double>(i) + 0.5) /
                                           static_cast<double>(n));
        } else {
            p = num_traits<T>::from_double(rng.uniform());
        }
        const Distrib<T> alt = sp.factory(dist_quantile(sp.param_dist, p));
        if (alt.disabled || alt.is_prior())
            throw InputError("prior_discretize: the factory returned an unusable Distribution");
        out.dists.push_back(alt);
        out.weights.push_back(wn);
    }
    return out;
}
 
/** `prior_discretize` with the defaults of `Prior.discretize`: 11 quadrature nodes. */
template <class T>
PriorDesign<T> prior_discretize(const Distrib<T>& d) {
    PriorRng rng(23000);
    return prior_discretize(d, kPriorDefaultNodes, "quadrature", rng);
}
 
/** E[X] = sum_i p_i E[X_i], MATLAB `Prior.getMean`. */
template <class T>
T prior_mean(const Distrib<T>& d) {
    const PriorDesign<T> g = prior_discretize(d);
    T m = num_traits<T>::from_int(0);
    for (std::size_t i = 0; i < g.dists.size(); ++i) m += T(g.weights[i] * g.dists[i].mean);
    return m;
}
 
/**
 * The SCV by the law of total variance, MATLAB `Prior.getSCV`.
 *
 * Var(X) = E[Var(X|D)] + Var(E[X|D]): the within-alternative variance plus the
 * spread of the alternative means. The second term is what makes a Prior's SCV
 * exceed the SCV of any of its alternatives.
 */
template <class T>
T prior_scv(const Distrib<T>& d) {
    const PriorDesign<T> g = prior_discretize(d);
    T e_mean = num_traits<T>::from_int(0), e_var = num_traits<T>::from_int(0),
      e_mean_sq = num_traits<T>::from_int(0);
    for (std::size_t i = 0; i < g.dists.size(); ++i) {
        const T m = g.dists[i].mean;
        const T v = T(g.dists[i].scv * m * m);
        e_mean += T(g.weights[i] * m);
        e_var += T(g.weights[i] * v);
        e_mean_sq += T(g.weights[i] * m * m);
    }
    if (!(e_mean > num_traits<T>::from_int(0)))
        throw NumericError("prior_scv: the prior-weighted mean is not positive");
    const T total_var = T(e_var + (e_mean_sq - e_mean * e_mean));
    return T(total_var / (e_mean * e_mean));
}
 
/**
 * The skewness of the mixture, MATLAB `Prior.getSkewness`.
 *
 * Each alternative's third central moment is shifted to the global mean by
 * E[(X_i - mu)^3] = m3_i + 3 v_i delta + delta^3, delta = m_i - mu, and the
 * shifted moments are averaged. Exact for the mixture, and it is the mixture
 * that `Distribution.getSkewness` is asked about.
 */
template <class T>
T prior_skewness(const Distrib<T>& d) {
    const T zero = num_traits<T>::from_int(0), three = num_traits<T>::from_int(3);
    const T mu = prior_mean(d);
    const T scv = prior_scv(d);
    const T sigma2 = T(scv * mu * mu);
    const double sigma = std::sqrt(num_traits<T>::to_double(sigma2));
    if (sigma < GlobalConstants::FineTol) return zero;
    const PriorDesign<T> g = prior_discretize(d);
    T third = zero;
    for (std::size_t i = 0; i < g.dists.size(); ++i) {
        const T mi = g.dists[i].mean;
        const T vi = T(g.dists[i].scv * mi * mi);
        const T m3i = T(dist_moment(g.dists[i], 3) - three * mi * vi - mi * mi * mi);
        const T delta = T(mi - mu);
        third += T(g.weights[i] * (m3i + three * vi * delta + delta * delta * delta));
    }
    return T(third / num_traits<T>::from_double(sigma * sigma * sigma));
}
 
/** F(t) = sum_i p_i F_i(t), MATLAB `Prior.evalCDF`. */
template <class T>
T prior_cdf(const Distrib<T>& d, const T& t) {
    const PriorDesign<T> g = prior_discretize(d);
    T f = num_traits<T>::from_int(0);
    for (std::size_t i = 0; i < g.dists.size(); ++i) f += T(g.weights[i] * dist_cdf(g.dists[i], t));
    return f;
}
 
/** L(s) = sum_i p_i L_i(s), MATLAB `Prior.evalLST`. */
template <class T>
T prior_lst(const Distrib<T>& d, const T& s) {
    const PriorDesign<T> g = prior_discretize(d);
    T l = num_traits<T>::from_int(0);
    for (std::size_t i = 0; i < g.dists.size(); ++i) l += T(g.weights[i] * dist_lst(g.dists[i], s));
    return l;
}
 
/**
 * Write the mixture moments onto a Prior, the counterpart of
 * `dist_refresh_moments` for the Markovian families.
 *
 * `prior_discrete` and `prior_continuous` leave `mean` and `scv` at their
 * defaults, because computing them costs a discretization of the parameter
 * density; the builders call this so that a struct dumped before SolverUQ has
 * run reports the epistemic mean rather than a zero.
 */
template <class T>
void prior_refresh_moments(Distrib<T>& d) {
    if (!d.is_prior()) return;
    d.mean = prior_mean(d);
    d.scv = prior_scv(d);
}
 
}  // namespace lang
}  // namespace line
 
#endif  // LINE_LANG_PRIOR_H