namespace gpr::laplace

Overview

namespace laplace {
 
// structs
 
struct PsisLaplaceOptions;
struct PsisLaplaceState;
struct SliceSamplerOptions;
struct SliceSamplerState;
struct SteinVariationalConfig;
 
// classes
 
class PsisLaplace;
class ScopedKernelState;
class SliceSampler;
class SteinVariational;
class TCovSecondOrder;
 
// global functions
 
double logGaussianDensity(const Eigen::VectorXd& delta, const Eigen::MatrixXd& cov_chol);
std::vector<double> selfNormalisedWeights(const std::vector<double>& log_w);
double estimateParetoKHat(const std::vector<double>& weights_descending, int M_pareto);
 
} // namespace laplace

Detailed Documentation

Global Functions

double logGaussianDensity(const Eigen::VectorXd& delta, const Eigen::MatrixXd& cov_chol)

Log-density of N(0, Sigma) at z, given the Cholesky factor L of Sigma (Sigma = L L^T) and the offset (theta - theta_map). Returns -0.5 (theta - theta_map)^T Sigma^{-1} (theta - theta_map)

• 0.5 log det Sigma - (p/2) log(2 pi). Computed via two triangular solves on L for the quadratic term and 2 * sum log diag(L) for the log-det.

std::vector<double> selfNormalisedWeights(const std::vector<double>& log_w)

Self-normalised importance weights. Numerically stable: subtracts max(log_w) before exp to avoid overflow. Returns weights summing to 1. If log_w is empty, returns empty.

double estimateParetoKHat(const std::vector<double>& weights_descending, int M_pareto)

Hill-style estimator of the Pareto shape k_hat on the upper-tail of the weight distribution. weights_descending must be sorted high to low (caller responsibility). M_pareto is the number of upper-tail samples to use; Vehtari 2024 recommends M_pareto = min(K/5, 3*sqrt(K)) but for K = 16 that is 3, which is too few for the production code path and the diagnostic; we use min(K-1, max(K/5, 8)) to keep the diagnostic honest. Returns 0 if M_pareto < 2 or all weights are equal.