Measurement with Uncertainty

Try it live on Compiler Explorer

Overview

Scientific and engineering measurements always carry uncertainty. This example demonstrates how to create a custom representation type that pairs measured values with their uncertainties and automatically propagates errors through calculations using the mp-units library.

The implementation provides first-order uncertainty propagation following standard error analysis formulas, making it suitable for experimental physics, metrology, and engineering applications. It's particularly relevant for systems like IAU (International Astronomical Union) where constants are defined with explicit uncertainties.

The example also demonstrates the Faster-than-lightspeed constants feature, where mathematical constants like π are treated as exact symbolic values that don't contribute to numerical uncertainty, ensuring that only measurement errors propagate through geometric calculations.

The measurement Class

Header File

The complete implementation is provided in a separate header file that can be reused across projects:

#pragma once

#ifdef MP_UNITS_IMPORT_STD
import std;
#else
#include <cmath>
#include <compare>  // IWYU pragma: export
#include <ostream>
#include <utility>
#endif

/**
 * @brief A representation type for physical measurements with uncertainties.
 *
 * This class represents a measured value with its associated uncertainty (standard deviation),
 * providing automatic uncertainty propagation through mathematical operations.
 *
 * @tparam T The underlying numeric type (e.g., double, float)
 *
 * **Uncertainty Propagation:**
 *
 * The class implements first-order uncertainty propagation (linear approximation) using
 * standard formulas from error analysis:
 *
 * - Addition/Subtraction: σ² = σ_x² + σ_y² (quadrature sum)
 * - Multiplication/Division: (σ/f)² = (σ_x/x)² + (σ_y/y)² (relative uncertainties in quadrature)
 * - Functions: σ_f = |df/dx| × σ_x (derivative times uncertainty)
 *
 * **Important Assumptions and Limitations:**
 *
 * 1. **Independent Variables:** All measurements are assumed to be statistically independent.
 *    Operations like `x - x` will give non-zero uncertainty (should be zero for perfectly
 *    correlated variables). For correlated measurements, a more sophisticated approach with
 *    covariance matrices would be needed.
 *
 * 2. **First-Order Approximation:** Uses linear approximation (first derivative only).
 *    This is accurate when uncertainties are small relative to the values. For large
 *    relative uncertainties or highly non-linear functions, higher-order terms may be needed.
 *
 * 3. **Gaussian Distribution:** Assumes uncertainties represent one standard deviation (σ)
 *    of normally distributed errors. Not suitable for systematic errors or non-Gaussian
 *    distributions without additional consideration.
 *
 * 4. **No Correlation Tracking:** The class does not track which measurements are derived
 *    from common sources, so it cannot detect or handle correlations automatically.
 *
 * **When to Use:**
 * - Combining independent measurements from different instruments
 * - Propagating random uncertainties through calculations
 * - Educational purposes and demonstrations
 *
 * **When NOT to Use:**
 * - When measurements are correlated (e.g., `f(x,x)` where x appears multiple times)
 * - When systematic uncertainties dominate
 * - When uncertainties are large relative to values (>10%)
 * - For Monte Carlo simulations (use direct sampling instead)
 *
 * @note This implementation is suitable for typical experimental physics and engineering
 *       calculations where independent measurements with small relative uncertainties are
 *       combined. For IAU astronomical constants and similar systems defined in terms of
 *       uncertainties, this provides appropriate uncertainty propagation.
 *
 * **Example:**
 * @code
 * auto length = measurement{10.0, 0.1} * m;  // 10.0 ± 0.1 m
 * auto width = measurement{5.0, 0.05} * m;   // 5.0 ± 0.05 m
 * auto area = length * width;                // 50.0 ± 0.71 m²
 * @endcode
 */
template<typename T>
class measurement {
public:
  using value_type = T;

  measurement() = default;

  /**
   * @brief Constructs a measurement with a value and uncertainty.
   *
   * @param val The measured value
   * @param err The uncertainty (absolute, will be converted to positive value)
   *
   * @note The uncertainty is stored as absolute value |err|, so negative uncertainties
   *       are automatically corrected.
   */
  // NOLINTNEXTLINE(bugprone-easily-swappable-parameters)
  constexpr explicit measurement(value_type val, const value_type& err = {}) :
      value_(std::move(val)), uncertainty_([&] {
        using namespace std;
        return abs(err);
      }())
  {
  }

  /// @brief Returns the central measured value
  [[nodiscard]] constexpr const value_type& value() const { return value_; }

  /// @brief Returns the absolute uncertainty (standard deviation)
  [[nodiscard]] constexpr const value_type& uncertainty() const { return uncertainty_; }

  /// @brief Returns the relative uncertainty (σ/x)
  [[nodiscard]] constexpr value_type relative_uncertainty() const { return uncertainty() / value(); }

  /// @brief Returns the lower bound of the uncertainty interval (value - σ)
  [[nodiscard]] constexpr value_type lower_bound() const { return value() - uncertainty(); }

  /// @brief Returns the upper bound of the uncertainty interval (value + σ)
  [[nodiscard]] constexpr value_type upper_bound() const { return value() + uncertainty(); }

  /// @brief Unary negation (uncertainty unchanged)
  [[nodiscard]] constexpr measurement operator-() const { return measurement(-value(), uncertainty()); }

  /**
   * @brief Addition of two measurements.
   * @note Assumes independent measurements. Formula: σ² = σ_x² + σ_y²
   */
  [[nodiscard]] friend constexpr measurement operator+(const measurement& lhs, const measurement& rhs)
  {
    using namespace std;
    return measurement(lhs.value() + rhs.value(), hypot(lhs.uncertainty(), rhs.uncertainty()));
  }

  /**
   * @brief Subtraction of two measurements.
   * @note Assumes independent measurements. Formula: σ² = σ_x² + σ_y²
   * @warning For correlated measurements (e.g., x - x), this will incorrectly give non-zero uncertainty.
   */
  [[nodiscard]] friend constexpr measurement operator-(const measurement& lhs, const measurement& rhs)
  {
    using namespace std;
    return measurement(lhs.value() - rhs.value(), hypot(lhs.uncertainty(), rhs.uncertainty()));
  }

  /**
   * @brief Multiplication of two measurements.
   * @note Assumes independent measurements. Formula: (σ_f/f)² = (σ_x/x)² + (σ_y/y)²
   */
  [[nodiscard]] friend constexpr measurement operator*(const measurement& lhs, const measurement& rhs)
  {
    const auto val = lhs.value() * rhs.value();
    using namespace std;
    return measurement(val, val * hypot(lhs.relative_uncertainty(), rhs.relative_uncertainty()));
  }

  /**
   * @brief Multiplication by an exact scalar.
   * @note The scalar is treated as exact (no uncertainty). Formula: σ_f = |k| × σ_x
   */
  [[nodiscard]] friend constexpr measurement operator*(const measurement& lhs,
                                                       const std::convertible_to<value_type> auto& value)
  {
    using namespace std;
    return measurement(static_cast<value_type>(lhs.value() * value),
                       static_cast<value_type>(abs(value) * lhs.uncertainty()));
  }

  /// @brief Multiplication by an exact scalar (commutative).
  [[nodiscard]] friend constexpr measurement operator*(const std::convertible_to<value_type> auto& value,
                                                       const measurement& rhs)
  {
    using namespace std;
    return measurement(value * rhs.value(), abs(value) * rhs.uncertainty());
  }

  /**
   * @brief Division of two measurements.
   * @note Assumes independent measurements. Formula: (σ_f/f)² = (σ_x/x)² + (σ_y/y)²
   */
  [[nodiscard]] friend constexpr measurement operator/(const measurement& lhs, const measurement& rhs)
  {
    const auto val = lhs.value() / rhs.value();
    using namespace std;
    return measurement(val, val * hypot(lhs.relative_uncertainty(), rhs.relative_uncertainty()));
  }

  /**
   * @brief Division by an exact scalar.
   * @note The scalar is treated as exact (no uncertainty). Formula: σ_f = σ_x / |k|
   */
  [[nodiscard]] friend constexpr measurement operator/(const measurement& lhs,
                                                       const std::convertible_to<value_type> auto& value)
  {
    using namespace std;
    return measurement(static_cast<value_type>(lhs.value()) / value,
                       static_cast<value_type>(lhs.uncertainty() / abs(value)));
  }

  /**
   * @brief Division of an exact scalar by a measurement.
   * @note Formula: σ_f = |f| × (σ_x / x)
   */
  [[nodiscard]] friend constexpr measurement operator/(const std::convertible_to<value_type> auto& value,
                                                       const measurement& rhs)
  {
    const auto val = static_cast<value_type>(value / rhs.value());
    using namespace std;
    return measurement(val, abs(val) * rhs.relative_uncertainty());
  }

  [[nodiscard]] constexpr auto operator<=>(const measurement&) const = default;

  friend std::ostream& operator<<(std::ostream& os, const measurement& v)
  {
    return os << v.value() << " ± " << v.uncertainty();
  }

  /**
   * @brief Absolute value of a measurement.
   * @note Uncertainty is preserved unchanged.
   */
  [[nodiscard]] friend constexpr measurement abs(const measurement& v)
    requires requires { abs(v.value()); } || requires { std::abs(v.value()); }
  {
    using std::abs;
    return measurement(abs(v.value()), v.uncertainty());
  }

  /**
   * @brief Power function with runtime exponent.
   *
   * @param base The measurement to raise to a power
   * @param exponent The exponent (exact value, no uncertainty)
   * @return measurement Result with propagated uncertainty
   *
   * @note Formula: if f = x^n, then σ_f = |n × x^(n-1) × σ_x| = |n × f/x × σ_x|
   * @note This is a first-order approximation valid for small relative uncertainties.
   */
  [[nodiscard]] friend constexpr measurement pow(const measurement& base, const value_type& exponent)
    requires requires { pow(base.value(), exponent); } || requires { std::pow(base.value(), exponent); }
  {
    using std::pow;
    using std::abs;
    const auto val = pow(base.value(), exponent);
    return measurement(val, abs(exponent * val / base.value() * base.uncertainty()));
  }

  /**
   * @brief Square root of a measurement.
   *
   * @note Formula: σ_f = σ_x / (2√x)
   * @note This is equivalent to pow(v, 0.5) but more efficient.
   */
  [[nodiscard]] friend constexpr measurement sqrt(const measurement& v)
    requires requires { sqrt(v.value()); } || requires { std::sqrt(v.value()); }
  {
    using std::sqrt;
    const auto val = sqrt(v.value());
    return measurement(val, v.uncertainty() / (value_type{2} * val));
  }

  /**
   * @brief Exponential function of a measurement.
   *
   * @note Formula: if f = exp(x), then σ_f = |exp(x) × σ_x| = |f × σ_x|
   * @warning Uncertainty grows exponentially with the value. For large x, this can lead
   *          to very large relative uncertainties where first-order approximation breaks down.
   */
  [[nodiscard]] friend constexpr measurement exp(const measurement& v)
    requires requires { exp(v.value()); } || requires { std::exp(v.value()); }
  {
    using std::exp;
    using std::abs;
    const auto val = exp(v.value());
    return measurement(val, abs(val * v.uncertainty()));
  }

  /**
   * @brief Natural logarithm of a measurement.
   *
   * @note Formula: if f = ln(x), then σ_f = |σ_x / x| = σ_x / x (for positive x)
   * @note Relative uncertainty in x becomes absolute uncertainty in ln(x).
   */
  [[nodiscard]] friend constexpr measurement log(const measurement& v)
    requires requires { log(v.value()); } || requires { std::log(v.value()); }
  {
    using std::log;
    using std::abs;
    return measurement(log(v.value()), abs(v.uncertainty() / v.value()));
  }

private:
  value_type value_{};
  value_type uncertainty_{};
};

Key implementation details

• Value and uncertainty storage: Stores both the measured value and its absolute uncertainty (σ)
• Automatic error propagation: All arithmetic operators implement proper uncertainty formulas
• Mathematical functions: Includes pow, sqrt, exp, log with correct uncertainty derivatives
• First-order approximation: Uses linear Taylor expansion for uncertainty propagation
• Independent variables assumption: Does not track correlations between measurements

Integration with mp-units

To use a custom type as a quantity representation, it must satisfy the RepresentationOf concept. The library verifies this at compile time:

measurement.cpp

static_assert(mp_units::RepresentationOf<measurement<double>, mp_units::quantity_character::real_scalar>);
static_assert(mp_units::RepresentationOf<measurement<double>, mp_units::quantity_character::vector>);

This allows measurement<T> to be used seamlessly with any mp-units quantity type.

Example Usage

Basic Operations

The example demonstrates various uncertainty propagation scenarios:

measurement.cpp

void example()
{
  using namespace mp_units;
  using namespace mp_units::si::unit_symbols;

  std::cout << "Mass of the Sun:        M_sun = " << measurement{19884, 2} * (mag_power<10, 26> * kg) << '\n';

  const auto acceleration = isq::acceleration(measurement{9.8, 0.1} * m / s2);
  const auto time = measurement{1.2, 0.1} * s;
  const QuantityOf<isq::velocity> auto velocity = acceleration * time;
  std::cout << "Velocity calculation:   V = " << acceleration << " * " << time << " = " << velocity << " = "
            << velocity.in(km / h) << '\n';

  const auto length = measurement{123., 1.} * m;
  std::cout << "Scalar multiplication:  d = 10 * " << length << " = " << 10 * length << '\n';

  const auto radius = measurement{5.0, 0.1} * m;
  const auto circumference = radius * (mag<2> * mag<π> * one);
  const auto area = pow<2>(radius) * (mag<π> * one);
  std::cout << "Radius:                 r = " << radius << '\n';
  std::cout << "Circular circumference: 2πr = " << circumference << " = " << circumference.in(m) << '\n';
  std::cout << "Circular area:          πr² = " << area << " = " << area.in(m2) << '\n';

  const auto area_measured = measurement{25.0, 1.0} * (mag<π> * m2);
  const auto radius_from_area = sqrt(area_measured / (mag<π> * one));
  std::cout << "Radius from area:       A = " << area_measured << " -> r = √(A/π) = " << radius_from_area << '\n';
}

This showcases:

1. Astronomical constants - Using magnitudes with uncertainties for large-scale values
2. Kinematic calculations - Multiplying acceleration by time with automatic unit conversions
3. Scalar multiplication - Exact scaling preserves relative uncertainty
4. Faster-than-lightspeed constants - Using mag<π> keeps π as an exact symbolic value, not a numeric approximation
5. Geometric calculations - Computing circumference and area where π contributes zero uncertainty
6. Inverse operations - Computing radius from area using sqrt, with π canceling out exactly

Sample Output:

Mass of the Sun:        M_sun = 19884 ± 2 (10²⁶ kg)
Velocity calculation:   V = 9.8 ± 0.1 m/s² * 1.2 ± 0.1 s = 11.76 ± 0.98732 m/s = 42.336 ± 3.55435 km/h
Scalar multiplication:  d = 10 * 123 ± 1 m = 1230 ± 10 m
Radius:                 r = 5 ± 0.1 m
Circular circumference: 2πr = 5 ± 0.1 (2 π) m = 31.4159 ± 0.628319 m
Circular area:          πr² = 25 ± 1 (π) m² = 78.5398 ± 3.14159 m²
Radius from area:       A = 25 ± 1 (π m²) -> r = √(A/π) = 5 ± 0.1 m

Uncertainty Propagation Analysis

Notice how uncertainties propagate through the calculations:

• Velocity: Combining two independent measurements (acceleration and time), both with ~1-2% relative uncertainty, yields a velocity with ~8.4% relative uncertainty (0.98732/11.76 ≈ 0.084)
• Scalar multiplication: Multiplying by exact scalar 10 preserves the relative uncertainty (~0.8%)
• Circumference (2πr): The factor 2π is exact (no uncertainty), so only the radius uncertainty propagates: σ_c = 2π × 0.1 m = 0.628... m
• Area (πr²): Squaring doubles the relative uncertainty (r: 2% → r²: 4%), while π remains exact: σ_A = π × |2r × σ_r| = π × 1 m² = 3.14... m²
• Round-trip verification: Computing radius from area with √(A/π) recovers the original 5 ± 0.1 m because π cancels exactly—this demonstrates the power of symbolic constants in avoiding accumulated numerical errors

Error Propagation Formulas

The implementation uses standard first-order uncertainty propagation:

Operation	Formula	Implementation
Addition: z = x + y	σ_z² = σ_x² + σ_y²	Quadrature sum (Pythagorean)
Subtraction: z = x - y	σ_z² = σ_x² + σ_y²	Same as addition
Multiplication: z = x × y	(σ_z/z)² = (σ_x/x)² + (σ_y/y)²	Relative uncertainties in quadrature
Division: z = x / y	(σ_z/z)² = (σ_x/x)² + (σ_y/y)²	Same as multiplication
Scalar multiply: z = k × x	σ_z = |k| × σ_x	Exact scaling
Power: z = x^n	σ_z = |n × z/x × σ_x|	Derivative method
Square root: z = √x	σ_z = σ_x / (2√x)	Special case of power
Exponential: z = exp(x)	σ_z = |z × σ_x|	Derivative method
Logarithm: z = ln(x)	σ_z = |σ_x / x|	Derivative method

Limitations and Caveats

Independent measurements assumption

The measurement class assumes all values are statistically independent. Operations like x - x will incorrectly give non-zero uncertainty. For correlated measurements, a more sophisticated approach with covariance matrices is needed.

First-order approximation

Uncertainty propagation uses linear approximation (first derivative only). This is accurate when relative uncertainties are small (<10%). For larger uncertainties or highly non-linear functions, higher-order terms may be significant.

When to use this implementation

• ✅ Combining independent measurements from different instruments
• ✅ Standard metrology and experimental physics calculations
• ✅ Systems like IAU where constants have defined uncertainties
• ✅ Small to moderate relative uncertainties (<10%)
• ❌ Correlated measurements (same source used multiple times)
• ❌ Large relative uncertainties (>10%)
• ❌ Systematic errors (requires different treatment)

Why This Matters

• Automatic Error Propagation: No manual uncertainty calculations needed—formulas are built into operators
• Type Safety: The dimensional analysis works with uncertainties seamlessly through mp-units
• Scientific Accuracy: Properly tracks measurement precision through complex calculations
• Standard Compliant: Follows ISO GUM (Guide to the Expression of Uncertainty in Measurement) principles
• Extensibility: Demonstrates how to integrate domain-specific numeric types with mp-units
• Real-World Applicability: Suitable for IAU astronomical constants, NIST physical constants, and laboratory measurements

This pattern is essential for scientific computing, metrology, laboratory measurements, and any application where measurement precision and traceability matter.