calibration_utils

Description of helper functions for calibration step implemented in the module bnn_for_14C_calibration.calibration_utils:

bnn_for_14C_calibration.calibration_utils

compute_HPD_regions(alpha, density=None, nb_intervals=1000, support_bounds=(0, 1), subdivision_components=None)

Compute the Highest Posterior Density (HPD) region for a univariate posterior distribution represented by a piecewise-constant density function on a regular subdivision.

Let \(f\) be a density of the posterior distribution. Then, the HPD region for credibility level (\(1 - \alpha\)) is the set: $$ \text{HPD} = \{ x : f(x) \ge k_{1-\alpha} \} $$ where \(k_{1-\alpha}\) is the smallest value such that the total posterior mass of the set \(\{ x : f(x) \ge k_{1-\alpha} \}\) is at least (\(1 - \alpha\)).

This function approximates the posterior density using midpoint densities on nb_intervals sub-intervals. The HPD region may be disconnected; the function returns all connected components of the HPD set.

Parameters:

Name	Type	Description	Default
alpha	float	Tail probability. The HPD region contains mass (1 - alpha). Must satisfy 0 <= alpha <= 1.	required
density	callable	A function density(dates: np.ndarray) -> np.ndarray computing the approximate posterior density at given points. Required if subdivision_components is None.	None
nb_intervals	int	Number of equal-length subintervals of the support on which the density is approximated as piecewise constant (based on midpoint evaluation). Required if subdivision_components is None.	1000
support_bounds	tuple of float	Tuple (min, max) specifying the bounds of the scaled date support. Only (0,1) is currently supported. Default is (0,1).	(0, 1)
subdivision_components	tuple or list of three numpy.ndarray	Precomputed components for the density approximation: - array of interval bounds (length N+1) - array of middle points (length N) - array of density values at middle points (length N) If provided, density and nb_intervals are ignored.	None

Returns:

Type	Description
dict	{ "calage_posterior_mode": float, "calage_posterior_mode_density": float, "connexe_HPD_intervals": list of [a,b] lists, "connexe_HPD_intervals_density": list of float, "HPD_threshold": float }

Raises:

Type	Description
NotImplementedError	If support_bounds is not (0, 1).
ValueError	If neither density nor subdivision_components are provided, or if subdivision_components does not contain exactly three arrays.

Notes

Piecewise-constant density approximation

When the user does not provide subdivision_components, the support is subdivided into nb_intervals equal-length intervals of width:
$$ h = \dfrac{1 - 0}{\text{nb_intervals}}. $$ The posterior density is then approximated as constant on each interval, equal to the value at its midpoint. The normalisation is computed via the sum of these midpoint densities, which is mathematically equivalent to a trapezoidal scheme specialised to constant-per-interval densities.

Computation of the HPD region

Let \((f_i)_{1 \le i \le N}\) be the midpoint densities and let \((I_i)_{1 \le i \le N}\) be their associated intervals. Sorting the \(f_i\) in decreasing order produces a sequence of density levels ordered from the most probable to the least probable regions of the posterior.

Define the scaled cumulative sum: $$ S_k = \dfrac{f_{(1)} + \cdots + f_{(k)}}{\sum_{i=1}^N f_i}, $$ where \((j)\) denotes the ordering from largest to smallest. The index \(k\) such that: $$ S_k \ge 1 - \alpha $$ determines the HPD density threshold: $$ k_{1-\alpha} = \dfrac{f_{(k)}}{\sum_{i=1}^N f_i}. $$ All intervals whose midpoint density is equal or greater than \(k_{1-\alpha}\) belong to the HPD region. Adjacent selected intervals are merged into connected components.

The function returns:
- the HPD threshold \(k_{1-\alpha}\),
- the connected HPD components expressed in terms of scaled dates (see minimax_scaling),
- the (scaled) posterior mode and its density.

Examples:

>>> f = lambda x: 2*(1 - x)
>>> res = compute_HPD_regions(alpha=0.1, density=f, nb_intervals=1000)
>>> "connexe_HPD_intervals" in res
True

Source code in src/bnn_for_14C_calibration/calibration_utils.py

def compute_HPD_regions(
    alpha: float,
    density: Optional[Callable[[np.ndarray], np.ndarray]] = None,
    nb_intervals: int = 1000,
    support_bounds: Tuple[float, float] = (0, 1),
    subdivision_components: Optional[
        Union[
            Tuple[np.ndarray, np.ndarray, np.ndarray],
            List[np.ndarray]
        ]
    ] = None
) -> Dict[str, Any]:
    """
    Compute the Highest Posterior Density (HPD) region for a univariate posterior
    distribution represented by a piecewise-constant density function 
    on a regular subdivision.

    Let $f$ be a density of the posterior distribution. 
    Then, the HPD region for credibility level ($1 - \\alpha$) is the set: 
    $$ 
        \\text{HPD} = \\\{ x : f(x) \ge k_{1-\\alpha} \\\}
    $$
    where $k_{1-\\alpha}$ is the smallest value such that the total posterior mass of the
    set $\\{ x : f(x) \ge k_{1-\\alpha} \\}$ is at least ($1 - \\alpha$).

    This function approximates the posterior density using midpoint densities on
    `nb_intervals` sub-intervals. The HPD region may be disconnected; the function returns
    all connected components of the HPD set.

    Parameters
    ----------
    alpha : float
        Tail probability. The HPD region contains mass (1 - alpha). Must satisfy 0 <= alpha <= 1.
    density : callable, optional
        A function `density(dates: np.ndarray) -> np.ndarray` computing the 
        approximate posterior density at given points.  
        Required if `subdivision_components` is None.
    nb_intervals : int, default=1000
        Number of equal-length subintervals of the support on which the density is
        approximated as piecewise constant (based on midpoint evaluation).  
        Required if `subdivision_components` is None.
    support_bounds : tuple of float, optional
        Tuple `(min, max)` specifying the bounds of the scaled date support. Only `(0,1)` is
        currently supported. Default is `(0,1)`.
    subdivision_components : tuple or list of three numpy.ndarray, optional
        Precomputed components for the density approximation:  
            - array of interval bounds (length N+1)  
            - array of middle points  (length N)  
            - array of density values at middle points  (length N)  
        If provided, `density` and `nb_intervals` are ignored.

    Returns
    -------
    dict
        {  
            "calage_posterior_mode": float,  
            "calage_posterior_mode_density": float,  
            "connexe_HPD_intervals": list of [a,b] lists,  
            "connexe_HPD_intervals_density": list of float,  
            "HPD_threshold": float  
        }

    Raises
    ------
    NotImplementedError
        If `support_bounds` is not `(0, 1)`.
    ValueError
        If neither `density` nor `subdivision_components` are provided, or if
        `subdivision_components` does not contain exactly three arrays.

    Notes
    -----
    **Piecewise-constant density approximation**

    When the user does not provide `subdivision_components`, the support is subdivided
    into `nb_intervals` equal-length intervals of width:  
    $$
        h = \dfrac{1 - 0}{\\text{nb_intervals}}.
    $$
    The posterior density is then approximated as *constant on each interval*, equal to
    the value at its midpoint. The normalisation is computed via the sum of these
    midpoint densities, which is mathematically equivalent to a trapezoidal scheme
    specialised to constant-per-interval densities.

    **Computation of the HPD region**

    Let $(f_i)_{1 \le i \le N}$ be the midpoint densities and let $(I_i)_{1 \le i \le N}$ be 
    their associated intervals.
    Sorting the $f_i$ in decreasing order produces a sequence of density levels ordered
    from the most probable to the least probable regions of the posterior.

    Define the scaled cumulative sum:
    $$
        S_k = \\dfrac{f_{(1)} + \cdots + f_{(k)}}{\sum_{i=1}^N f_i},
    $$
    where $(j)$ denotes the ordering from largest to smallest. The index $k$ such that:
    $$
        S_k \ge 1 - \\alpha
    $$
    determines the HPD density threshold:
    $$
        k_{1-\\alpha} = \dfrac{f_{(k)}}{\sum_{i=1}^N f_i}.
    $$
    All intervals whose midpoint density is equal or greater than $k_{1-\\alpha}$ belong to 
    the HPD region. Adjacent selected intervals are merged into connected components.

    The function returns:  
        - the HPD threshold $k_{1-\\alpha}$,  
        - the connected HPD components expressed in terms of scaled dates
            (see `minimax_scaling`),  
        - the (scaled) posterior mode and its density.  

    Examples
    --------
    >>> f = lambda x: 2*(1 - x)
    >>> res = compute_HPD_regions(alpha=0.1, density=f, nb_intervals=1000)
    >>> "connexe_HPD_intervals" in res
    True
    """

    # traitement du support et contrôle des arguments fournis
    if support_bounds != (0,1) :
        raise NotImplementedError(
            "Only the default support (0,1) is currently supported"
        )

    if density == None and subdivision_components == None :
        raise ValueError(
            "At least one of 'density' or 'subdivision_components' must be provided (not None)"
        )

    if subdivision_components != None and len(subdivision_components) != 3 :
        raise ValueError(
            "'subdivision_components' must be a tuple or list of length 3 containing arrays: interval bounds, middle points, and densities at middle points"
        )

    if subdivision_components != None : 

        intervals_bounds = subdivision_components[0]
        middle_points = subdivision_components[1]
        middle_points_density = subdivision_components[2]
        # nb_intervals = len(middle_points)

    else :

        # subdivision du support en nb_intervals : calcul des bornes des sous-intervalles
        intervals_bounds = np.linspace(support_bounds[0], support_bounds[1], nb_intervals+1, dtype=np.float64)

        # évaluation de la densité aux points milieu
        # midle_points = (support_bounds[1] - support_bounds[0])/(2*nb_intervals) + intervals_bounds[:-1]
        middle_points = (intervals_bounds[:-1] + intervals_bounds[1:])/2
        middle_points_density = density(middle_points)

    # on reordonne les intervalles (les densités) suivant la valeur de leurs densités, dans l'ordre décroissante
    sorted_desc_index = np.argsort(middle_points_density)[::-1]
    sorted_desc_middle_points_density = middle_points_density[sorted_desc_index]

    # on calcule la densité cumulée des densités ainsi reordonnéees et on les renormalise pour 
    # pouvoir avoir 1 comme dernier élément comme dans unvecteur de fonction de répartition
    scaling_weight = sorted_desc_middle_points_density.sum()
    sorted_desc_middle_points_density_cumsum_scaled = sorted_desc_middle_points_density.cumsum()/scaling_weight

    # au cas où il y a des valeurs > 1 dans sorted_desc_middle_points_density_cumsum_scaled (à cause des erreurs d'arrondis), on les met à 1
    sorted_desc_middle_points_density_cumsum_scaled = np.where(sorted_desc_middle_points_density_cumsum_scaled > 1., 1., sorted_desc_middle_points_density_cumsum_scaled)

    # on calcule le mode a posteriori
    # (TODO : voir plus tard si nécessaire de calculer plusieurs modes de même densité HPD éventuellement)
    calage_posterior_mode_density = middle_points_density[sorted_desc_index[0]]/scaling_weight
    calage_posterior_mode = middle_points[sorted_desc_index[0]]


    # on peut alors regarder à partir de quand on dépasse le seuil de 1 - alpha
    where_idx = np.where(sorted_desc_middle_points_density_cumsum_scaled >= 1 - alpha)[0][0]
    k_1_alpha = sorted_desc_middle_points_density[where_idx]/scaling_weight # = quantile d'ordre alpha du vecteur des densités

    # on recupère alors les (indices des) intervalles formant la région HPD sélectionnés parmi les nb_intervals : 
    # les indices sont reordonnés du plus petit au plus grand
    selected_intervals_idx = np.sort(sorted_desc_index[:(where_idx+1)])

    # maintenant, il suffit de déterminer les intervalles HPD sélectionnés, tout en regroupant les parties connexes 
    # et en calculant la densité de chaque partie connexe ainsi constituée
    connexe_intervals = []
    connexe_intervals_density = []
    l = len(selected_intervals_idx)
    i = 0
    while i < l :
        first = selected_intervals_idx[i]
        last = first
        j = i+1
        while j < l :
            current = selected_intervals_idx[j]
            if current - last == 1 :
                last = current
                j = j+1
            else :
                break
        if first == last :
            connexe_intervals.append([intervals_bounds[first], intervals_bounds[first+1]])
            connexe_intervals_density.append(middle_points_density[first]/scaling_weight)
        else :
            connexe_intervals.append([intervals_bounds[first], intervals_bounds[last+1]])
            connexe_intervals_density.append(middle_points_density[first:(last+1)].sum()/scaling_weight)
        i = j

    # on peut alors retourner la région HPD sous formes des parties connexes ainsi que les densités
    # de différentes parties connexes
    # en bonus, on ajoute le seuil k_1_alpha permettant de déterminer la région HPD
    # on ajoute aussi le mode calage_mode le plus plaussible : la date calibrée la plus probable a posteriori
    return {
            "calage_posterior_mode" : calage_posterior_mode,
            "calage_posterior_mode_density" : calage_posterior_mode_density,
            "connexe_HPD_intervals" : connexe_intervals, 
            "connexe_HPD_intervals_density" : connexe_intervals_density, 
            "HPD_threshold" : k_1_alpha
        }

mono_cal_date_approx_cumulative_fct(density=None, nb_intervals=1000, support_bounds=(0, 1), subdivision_components=None)

Approximate the posterior cumulative distribution function (CDF) for a single calibrated radiocarbon date.

This function builds a continuous CDF from either:
- a callable density function density, or
- precomputed subdivision components (interval bounds, middle points, densities at middle points).

The continuity of the resulting CDF is a direct mathematical consequence of the sampling strategy used for generating posterior samples: each interval is chosen with probability proportional to its middle-point density, and points within the interval are drawn uniformly. This results in a CDF that increases linearly within each interval.

Parameters:

Name	Type	Description	Default
density	callable	A function density(dates: np.ndarray) -> np.ndarray computing the approximate posterior density at given points. Required if subdivision_components is None.	None
nb_intervals	int	Number of subintervals to discretize the support for the density approximation. Overrided internally if subdivision_components is provided.	1000
support_bounds	tuple of float	Tuple (min, max) specifying the bounds of the scaled date support. Only (0,1) is currently supported. Default is (0,1).	(0, 1)
subdivision_components	tuple or list of three numpy.ndarray	Precomputed components for the density approximation: - array of interval bounds (length N+1) - array of middle points (length N) - array of density values at middle points (length N) If provided, density is ignored and nb_intervals is overrided.	None

Returns:

Name	Type	Description
cumulative_density	callable	A function CDF(date: float) -> float returning the approximate posterior cumulative probability for the given date.

Raises:

Type	Description
NotImplementedError	If support_bounds is not (0, 1).
ValueError	If neither density nor subdivision_components are provided, or if subdivision_components does not contain exactly three arrays.

Notes

• Let the support \([a,b]\) be subdivided into N intervals \([x_j, x_{j+1}]\) with middle points \(m_j = \frac{x_j + x_{j+1}}{2}\), and let \(f_j = density(m_j)\). The approximate \(CDF\) at a point \(d \in [x_j, x_{j+1}]\) is: $$ CDF(d) = \sum_{i=1}^{j-1} \dfrac{f_i}{\sum_{k=1}^N f_k} + \dfrac{f_j}{\sum_{k=1}^N f_k} \dfrac{d - x_j}{x_{j+1} - x_j} $$
  where the first term sums contributions from previous intervals, and the second term accounts for the uniform distribution inside the current interval.
• Continuity of the CDF arises naturally from the uniform distribution inside intervals.
• This approach provides a simple and fast approximation suitable for sampling posterior dates. For N → ∞, the discrete sum converges to the integral of the continuous piecewise density function.

Examples:

>>> density_fn = mono_cal_date_approx_density(
...     mesure=0.954,
...     lab_error=0.002,
...     bnn_model=my_trained_bnn,
...     nb_curves=200
... )
>>> cdf_fn = mono_cal_date_approx_cumulative_fct(density=density_fn, nb_intervals=1000)
>>> scaled_date = 0.25
>>> cdf_value = cdf_fn(scaled_date)
>>> isinstance(cdf_value, float)
True

Source code in src/bnn_for_14C_calibration/calibration_utils.py

def mono_cal_date_approx_cumulative_fct(
    density: Optional[Callable[[np.ndarray], np.ndarray]] = None,
    nb_intervals: int = 1000,
    support_bounds: Tuple[float, float] = (0, 1),
    subdivision_components: Optional[Union[Tuple[np.ndarray, np.ndarray, np.ndarray], List[np.ndarray]]] = None
) -> Callable[[float], float]:
    """
    Approximate the posterior cumulative distribution function (CDF) for a single calibrated radiocarbon date.

    This function builds a continuous CDF from either:  
      - a callable density function `density`, or  
      - precomputed subdivision components (interval bounds, middle points, densities at middle points).

    The continuity of the resulting CDF is a direct mathematical consequence of the
    sampling strategy used for generating posterior samples: each interval is chosen
    with probability proportional to its middle-point density, and points within the
    interval are drawn uniformly. This results in a CDF that increases linearly
    within each interval.

    Parameters
    ----------
    density : callable, optional
        A function `density(dates: np.ndarray) -> np.ndarray` computing the 
        approximate posterior density at given points. Required if `subdivision_components` is None.
    nb_intervals : int, default=1000
        Number of subintervals to discretize the support for the density approximation.
        Overrided internally if `subdivision_components` is provided.
    support_bounds : tuple of float, optional
        Tuple `(min, max)` specifying the bounds of the scaled date support. Only `(0,1)` is
        currently supported. Default is `(0,1)`.
    subdivision_components : tuple or list of three numpy.ndarray, optional
        Precomputed components for the density approximation:  
            - array of interval bounds (length N+1)  
            - array of middle points  (length N)  
            - array of density values at middle points  (length N)  
        If provided, `density` is ignored and `nb_intervals` is overrided.

    Returns
    -------
    cumulative_density : callable
        A function `CDF(date: float) -> float` returning the approximate
        posterior cumulative probability for the given date.

    Raises
    ------
    NotImplementedError
        If `support_bounds` is not `(0, 1)`.
    ValueError
        If neither `density` nor `subdivision_components` are provided, or if
        `subdivision_components` does not contain exactly three arrays.

    Notes
    -----
    - Let the support $[a,b]$ be subdivided into N intervals $[x_j, x_{j+1}]$ with
      middle points $m_j = \\frac{x_j + x_{j+1}}{2}$, and let $f_j = density(m_j)$.
      The approximate $CDF$ at a point $d \in [x_j, x_{j+1}]$ is:
        $$
        CDF(d) = \sum_{i=1}^{j-1} \\dfrac{f_i}{\sum_{k=1}^N f_k}  + 
        \\dfrac{f_j}{\sum_{k=1}^N f_k} \\dfrac{d - x_j}{x_{j+1} - x_j}
        $$  
      where the first term sums contributions from previous intervals, and the second
      term accounts for the uniform distribution inside the current interval.  
    - Continuity of the CDF arises naturally from the uniform distribution inside
      intervals.  
    - This approach provides a simple and fast approximation suitable for sampling
      posterior dates. For N → ∞, the discrete sum converges to the integral
      of the continuous piecewise density function.

    Examples
    --------
    >>> density_fn = mono_cal_date_approx_density(
    ...     mesure=0.954,
    ...     lab_error=0.002,
    ...     bnn_model=my_trained_bnn,
    ...     nb_curves=200
    ... )
    >>> cdf_fn = mono_cal_date_approx_cumulative_fct(density=density_fn, nb_intervals=1000)
    >>> scaled_date = 0.25
    >>> cdf_value = cdf_fn(scaled_date)
    >>> isinstance(cdf_value, float)
    True

    """

    # traitement du support et contrôle des arguments fournis
    if support_bounds != (0,1) :
        raise NotImplementedError(
            "Only the default support (0,1) is currently supported"
        )

    if density == None and subdivision_components == None :
        raise ValueError(
            "At least one of 'density' or 'subdivision_components' must be provided (not None)"
        )

    if subdivision_components != None and len(subdivision_components) != 3 :
        raise ValueError(
            "'subdivision_components' must be a tuple or list of length 3 containing arrays: interval bounds, middle points, and densities at middle points"
        )

    if subdivision_components != None : 

        intervals_bounds = subdivision_components[0]
        middle_points = subdivision_components[1]
        middle_points_density = subdivision_components[2]
        nb_intervals = len(middle_points)

    else :

        # subdivision du support en nb_intervals : calcul des bornes des sous-intervalles
        intervals_bounds = np.linspace(support_bounds[0], support_bounds[1], nb_intervals+1, dtype=np.float64)

        # évaluation de la densité aux points milieu
        # midle_points = (support_bounds[1] - support_bounds[0])/(2*nb_intervals) + intervals_bounds[:-1]
        middle_points = (intervals_bounds[:-1] + intervals_bounds[1:])/2
        middle_points_density = density(middle_points)

    # la borne inférieure de l'intervalle qui contient la date d sur laquelle évaluer la fonction
    lower_test = lambda d : intervals_bounds[:-1] <= d 
    idx = lambda d : np.where(lower_test(d))[0][-1] #-1 car c'est la dernière borne pour laquelle le test ci-dessus vaut true
    # NB : np.where renvoie un tuple d'array (ici c'est un tuple avec un seul array) et c'est le premier élément du tuple qui nous intéresse ici

    # calcul de la densité cumulée au point d 
    h = (support_bounds[1] - support_bounds[0])/nb_intervals
    cumulative_density = lambda d : (middle_points_density[idx(d)]*(d - intervals_bounds[idx(d)])/h + middle_points_density[:idx(d)].sum())/middle_points_density.sum()

    return cumulative_density

mono_cal_date_approx_density(mesure, lab_error, bnn_model, nb_curves=100, prior_density='default', batch_size=None)

Approximate the posterior density function for a single radiocarbon date calibration.

This function computes an approximate posterior density for a given measured radiocarbon date using a trained Bayesian Neural Network (BNN) model.
The approximation is done by sampling multiple stochastic realizations of the BNN predictions and combining them with a likelihood term based on the laboratory measurement error.

Parameters:

Name	Type	Description	Default
mesure	float	The measured radiocarbon age (expressed in the F\(^{14}\)C domain).	required
lab_error	float	The measurement uncertainty (standard deviation) associated with the lab measurement (also in the F\(^{14}\)C domain).	required
bnn_model	object	The trained Bayesian Neural Network model used to estimate the predictive distribution. Must be compatible with bnn_make_predictions_. This model must be an estimate of the radiocarbon calibration curve in the F\(^{14}\)C domain.	required
nb_curves	int	Number of stochastic realizations (Monte Carlo samples) to use for approximating the BNN predictive distribution. Default is 100.	100
prior_density	(default, callable)	Prior probability density over the calendar dates' domain. - "default": a uniform prior over [0, 1], the range of the scaled calendar dates. - callable: a custom prior density function of the form f(dates) → np.ndarray. Default is "default" (the only possibility supported presently).	"default"
batch_size	int	Batch size for model predictions, passed to the internal bnn_make_predictions_ function. Default is None.	None

Returns:

Name	Type	Description
density	callable	A function density(dates: np.ndarray) -> np.ndarray that computes the (unnormalized) approximate posterior density of calibrated dates. The dates given to this function are first scaled using minimax_scaling function. The density is known up to a normalization constant.

Notes

• The posterior density is proportional to the product of the prior and the likelihood:
  \( p(d|m) ∝ p(d) × E_{BNN}[ \exp(-(m - F^{14}C(d))^2 / (2σ^2)) ] \).
  This expectation is approximated by averaging over multiple stochastic predictions of the BNN model.
• The uniform prior currently assumes that the scaled ages' calibration domain is [0, 1];
  future implementations should allow to replace this with the actual domain limits or use other values to approximate these limits (e.g. training data domain bounds).
  Another implementation improvement may be to make it possible to handle the use of a callable which gives a custom prior density function of the form f(dates) → np.ndarray.
• The output density is not normalized; normalization must be handled externally if necessary (e.g., via numerical integration).

Examples:

>>> density_fn = mono_cal_date_approx_density(
...     mesure=0.954,
...     lab_error=0.002,
...     bnn_model=my_trained_bnn,
...     nb_curves=200
... )
>>> ages = np.linspace(0, 1, 500)
>>> posterior_vals = density_fn(ages)
>>> posterior_vals.shape
(500,)

Source code in src/bnn_for_14C_calibration/calibration_utils.py

def mono_cal_date_approx_density(
    mesure: float,
    lab_error: float,
    bnn_model: object,
    nb_curves: int = 100,
    prior_density: Union[str, Callable[[np.ndarray], np.ndarray]] = "default",
    batch_size: Optional[int] = None
) -> Callable[[np.ndarray], np.ndarray]:
    """
    Approximate the posterior density function for a single radiocarbon date calibration.

    This function computes an approximate posterior density for a given measured 
    radiocarbon date using a trained Bayesian Neural Network (BNN) model.  
    The approximation is done by sampling multiple stochastic realizations of 
    the BNN predictions and combining them with a likelihood term based on 
    the laboratory measurement error.

    Parameters
    ----------
    mesure : float
        The measured radiocarbon age (expressed in the F$^{14}$C domain).
    lab_error : float
        The measurement uncertainty (standard deviation) associated with the lab measurement 
        (also in the F$^{14}$C domain).
    bnn_model : object
        The trained Bayesian Neural Network model used to estimate the predictive distribution.  
        Must be compatible with `bnn_make_predictions_`.  
        This model must be an estimate of the radiocarbon calibration curve in the F$^{14}$C domain.
    nb_curves : int, optional
        Number of stochastic realizations (Monte Carlo samples) to use for 
        approximating the BNN predictive distribution.  
        Default is `100`.
    prior_density : {"default", callable}, optional
        Prior probability density over the calendar dates' domain.  
        - `"default"`: a uniform prior over `[0, 1]`, the range of the scaled calendar dates.  
        - `callable`: a custom prior density function of the form `f(dates) → np.ndarray`.  
        Default is `"default"` (the only possibility supported presently).
    batch_size : int, optional
        Batch size for model predictions, passed to the internal 
        `bnn_make_predictions_` function.  
        Default is `None`.

    Returns
    -------
    density : callable
        A function `density(dates: np.ndarray) -> np.ndarray` that computes 
        the (unnormalized) approximate posterior density of calibrated dates.  
        The dates given to this function are first scaled using `minimax_scaling` function.  
        The density is known up to a normalization constant.

    Notes
    -----
    - The posterior density is proportional to the product of the prior and the likelihood:  
      \\( p(d|m) ∝ p(d) × E_{BNN}[ \\exp(-(m - F^{14}C(d))^2 / (2σ^2)) ] \\).  
      This expectation is approximated by averaging over multiple stochastic 
      predictions of the BNN model.
    - The uniform prior currently assumes that the scaled ages' calibration domain is `[0, 1]`;  
      future implementations should allow to replace this with the actual domain limits or use
      other values to approximate these limits (e.g. training data domain bounds).  
      Another implementation improvement may be to make it possible to handle the use of a 
      `callable` which gives a custom prior density function of the form `f(dates) → np.ndarray`.  
    - The output density is **not normalized**; normalization must be handled externally 
      if necessary (e.g., via numerical integration).

    Examples
    --------
    >>> density_fn = mono_cal_date_approx_density(
    ...     mesure=0.954,
    ...     lab_error=0.002,
    ...     bnn_model=my_trained_bnn,
    ...     nb_curves=200
    ... )
    >>> ages = np.linspace(0, 1, 500)
    >>> posterior_vals = density_fn(ages)
    >>> posterior_vals.shape
    (500,)

    """

    # traitement de la densité à priori :
    if prior_density == "default":
        # à remplacer par min_Xtrain ou min_Xtrain_val ou min_Xtrain_val_test plus tard suivant le cas
        # ou date minimale globale possible pour la calibration
        support_lower_bound = 0.0
        # à remplacer par max_Xtrain ou max_Xtrain_val ou max_Xtrain_val_test plus tard suivant le cas
        # ou date maximale globale possible pour la calibration
        support_upper_bound = 1.0

        prior_density = lambda d: np.float64(
            (support_lower_bound <= d) * (d <= support_upper_bound)
        ) / (support_upper_bound - support_lower_bound)

    else:
        raise NotImplementedError(
            "Custom prior densities are not yet supported."
        )

    # predictions avec le modèle
    predicted = lambda d: bnn_make_predictions_(
        bnn_model=bnn_model,
        X_test=d.reshape((-1, 1)),
        iterations=nb_curves,
        batch_size=batch_size,
    )

    # densité approchée (connue à une constante près)
    density = lambda d: prior_density(d) * np.exp(
        -(mesure - predicted(d)) ** 2 / (2 * lab_error**2)
    ).mean(axis=1, dtype=np.float64) / (lab_error * np.sqrt(2 * np.pi))

    return density

mono_cal_date_approx_density_sample(density=None, nb_intervals=1000, support_bounds=(0.0, 1.0), subdivision_components=None, sample_size=1)

Draw samples from an univariate, unnormalized posterior density using a piecewise-constant approximation over a regular grid.

This function is typically used in Bayesian radiocarbon calibration to draw samples from the approximate posterior distribution of a single calibrated date, when the posterior density is only available up to a multiplicative constant or when direct numerical integration is not desirable.

Parameters:

Name	Type	Description	Default
density	callable or None	A function evaluating the unnormalized posterior density on an array of points. Required unless subdivision_components is supplied.	None
nb_intervals	int	Number of subintervals defining the grid approximation of the density.	1000
support_bounds	tuple of float	Lower and upper bounds of the support. Only (0, 1) is currently supported.	(0.0, 1.0)
subdivision_components	tuple of numpy.ndarray or None	Optional tuple (interval_bounds, midpoints, midpoint_densities). If provided, these arrays are reused directly.	None
sample_size	int	Number of posterior samples to generate.	1

Returns:

Type	Description
tuple	A tuple (d, unnorm_prob, norm_prob) with: - d : ndarray of shape (sample_size,), the generated samples. - unnorm_prob : ndarray of shape (sample_size,), the unnormalized density values at the selected midpoints. - norm_prob : ndarray of shape (sample_size,), the normalized discrete probabilities associated with the chosen intervals (see the notes below).

Notes

1. Posterior density availability

The function assumes that the posterior density `p(d | m)` is only known
through an unnormalized function:

    f(d) ∝ p(d | m)

This is the case when the density results from Monte Carlo averaging over a
Bayesian Neural Network (BNN), where:

    f(d) = p(d) × E_BNN[ exp(-(m - F¹⁴C(d))² / (2σ²)) ] / (σ √(2π))

Since the density is unnormalized, classical continuous inversion sampling 
is not possible. Instead, a piecewise-constant discretization is used.

2. Numerical approximation

The interval [0, 1] is subdivided into `nb_intervals` equal subintervals.
On each subinterval, the density is approximated by its midpoint value:

    f(d) ≈ f(d_j*)   for d in interval j

yielding discrete weights:

    p_j = f(d_j*) / Σ_k f(d_k*)

which define a categorical distribution over the intervals.

3. Sampling algorithm

Sampling is performed as follows:

1. Compute (or reuse) midpoints `d_j*` and their unnormalized densities.
2. Normalize these densities to obtain probabilities over subintervals.
3. Draw an interval index J according to these probabilities.
4. Draw a uniform sample on the chosen interval:

       d ~ Uniform(interval_bounds[J-1], interval_bounds[J])

This yields samples approximately distributed according to the target
posterior density.

4. Precomputed subdivision

If `subdivision_components = (bounds, midpoints, densities)` is supplied,
the function skips all density evaluations and directly reuses the
piecewise-constant representation.  
This is useful when repeatedly sampling from the same density is needed, 
e.g. inside an MCMC procedure.

5. Support and limitations

- Only the default support (0, 1) is currently implemented.
- The method is *univariate*.  
  Multivariate posterior sampling must instead rely on MCMC (e.g. 
  Metropolis–Hastings within Gibbs), since grid-based density approximation 
  in higher dimension is impractical due to the curse of dimensionality.
- The density is never normalized by continuous integration (on purpose),
  this yields a converging approximation of f(d) when the number of subintervals,
  `nb_intervals`, tends to infinity.

Examples:

>>> density = lambda x: np.exp(-(x-0.4)**2 / 0.01)
>>> d, f_unnorm, f_norm = mono_cal_date_approx_density_sample(
...     density=density, 
...     sample_size=5
... )

Source code in src/bnn_for_14C_calibration/calibration_utils.py

def mono_cal_date_approx_density_sample(
    density: Optional[Callable[[np.ndarray], np.ndarray]] = None,
    nb_intervals: int = 1000,
    support_bounds: Tuple[float, float] = (0.0, 1.0),
    subdivision_components: Optional[
        Tuple[np.ndarray, np.ndarray, np.ndarray]
    ] = None,
    sample_size: int = 1
) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    """
    Draw samples from an univariate, unnormalized posterior density using a 
    piecewise-constant approximation over a regular grid.

    This function is typically used in Bayesian radiocarbon calibration to draw
    samples from the approximate posterior distribution of a **single** calibrated
    date, when the posterior density is only available up to a multiplicative
    constant or when *direct* numerical integration is not desirable.

    Parameters
    ----------
    density : callable or None
        A function evaluating the unnormalized posterior density on an array
        of points. Required unless `subdivision_components` is supplied.

    nb_intervals : int
        Number of subintervals defining the grid approximation of the density.

    support_bounds : tuple of float
        Lower and upper bounds of the support. Only (0, 1) is currently supported.

    subdivision_components : tuple of numpy.ndarray or None
        Optional tuple (interval_bounds, midpoints, midpoint_densities).
        If provided, these arrays are reused directly.

    sample_size : int
        Number of posterior samples to generate.

    Returns
    -------
    tuple
        A tuple (d, unnorm_prob, norm_prob) with:  
            - d : `ndarray of shape (sample_size,)`, 
                the generated samples.  
            - unnorm_prob : `ndarray of shape (sample_size,)`, 
                the unnormalized density values at the selected midpoints.  
            - norm_prob : `ndarray of shape (sample_size,)`, 
                the normalized discrete probabilities associated with the chosen intervals 
                (see the notes below).  

    Notes
    -----
    **1. Posterior density availability**

        The function assumes that the posterior density `p(d | m)` is only known
        through an unnormalized function:

            f(d) ∝ p(d | m)

        This is the case when the density results from Monte Carlo averaging over a
        Bayesian Neural Network (BNN), where:

            f(d) = p(d) × E_BNN[ exp(-(m - F¹⁴C(d))² / (2σ²)) ] / (σ √(2π))

        Since the density is unnormalized, classical continuous inversion sampling 
        is not possible. Instead, a piecewise-constant discretization is used.

    **2. Numerical approximation**

        The interval [0, 1] is subdivided into `nb_intervals` equal subintervals.
        On each subinterval, the density is approximated by its midpoint value:

            f(d) ≈ f(d_j*)   for d in interval j

        yielding discrete weights:

            p_j = f(d_j*) / Σ_k f(d_k*)

        which define a categorical distribution over the intervals.

    **3. Sampling algorithm**

        Sampling is performed as follows:

        1. Compute (or reuse) midpoints `d_j*` and their unnormalized densities.
        2. Normalize these densities to obtain probabilities over subintervals.
        3. Draw an interval index J according to these probabilities.
        4. Draw a uniform sample on the chosen interval:

               d ~ Uniform(interval_bounds[J-1], interval_bounds[J])

        This yields samples approximately distributed according to the target
        posterior density.

    **4. Precomputed subdivision**  

        If `subdivision_components = (bounds, midpoints, densities)` is supplied,
        the function skips all density evaluations and directly reuses the
        piecewise-constant representation.  
        This is useful when repeatedly sampling from the same density is needed, 
        e.g. inside an MCMC procedure.

    **5. Support and limitations**

        - Only the default support (0, 1) is currently implemented.
        - The method is *univariate*.  
          Multivariate posterior sampling must instead rely on MCMC (e.g. 
          Metropolis–Hastings within Gibbs), since grid-based density approximation 
          in higher dimension is impractical due to the curse of dimensionality.
        - The density is never normalized by continuous integration (on purpose),
          this yields a converging approximation of f(d) when the number of subintervals,
          `nb_intervals`, tends to infinity.

    Examples
    --------
    >>> density = lambda x: np.exp(-(x-0.4)**2 / 0.01)
    >>> d, f_unnorm, f_norm = mono_cal_date_approx_density_sample(
    ...     density=density, 
    ...     sample_size=5
    ... )
    """

    # traitement du support et contrôle des arguments fournis
    if support_bounds != (0,1) :
        raise NotImplementedError(
            "Only support (0,1) is currently supported."
        )

    if density == None and subdivision_components == None :
        raise ValueError(
            "At least one of 'density' or 'subdivision_components' must be provided."
        )

    if subdivision_components != None and len(subdivision_components) != 3 :
        raise ValueError(
            "'subdivision_components' must be a tuple of three arrays: bounds, midpoints, densities."
        )

    if subdivision_components != None : 

        intervals_bounds = subdivision_components[0]
        middle_points = subdivision_components[1]
        middle_points_density = subdivision_components[2]
        nb_intervals = len(middle_points)

    else :

        # subdivision du support en nb_intervals : calcul des bornes des sous-intervalles
        intervals_bounds = np.linspace(support_bounds[0], support_bounds[1], nb_intervals+1, dtype=np.float64)

        # évaluation de la densité aux points milieu
        # midle_points = (support_bounds[1] - support_bounds[0])/(2*nb_intervals) + intervals_bounds[:-1]
        middle_points = (intervals_bounds[:-1] + intervals_bounds[1:])/2
        middle_points_density = density(middle_points)

    # pour tirer une date d suivant la densité voulue, on procède comme suit :
    #     *) tirer un indice j dans dans {1, 2, ..., nb_intervals} muni de la probabilité 
    #         middle_points_density/middle_points_density.sum() (on peut se contenter de la 
    #         probabilité non normalisée middle_points_density)
    #     *) tirer d suivant la loi uniforme sur le j ième intervalle [intervals_bounds[j-1], intervals_bounds[j]]
    rng = np.random.default_rng()
    probabilities = middle_points_density/middle_points_density.sum()
    # while probabilities.sum() != 1 :
    #     probabilities = probabilities/probabilities.sum()
    j = 1 + rng.choice(nb_intervals, size = sample_size, p = probabilities) # le +1 permet d'avoir j entre 1 et nb_intervals au lieu de 0 et nb_intervals - 1
    u = rng.random(size = sample_size) # loi uniforme sur [0,1]
    d = (intervals_bounds[j] - intervals_bounds[j-1]) * u + intervals_bounds[j-1] # équivalent aussi (support_bounds[1] - support_bounds[0])/nb_intervals * u + intervals_bounds[j-1]

    # on retourne d et sa probabilité (non normalisée et normalisée)
    return d, middle_points_density[j-1], probabilities[j-1]

mono_cal_date_approx_vect_cumulative_fct(density=None, nb_intervals=1000, support_bounds=(0, 1), subdivision_components=None)

Vectorized version of mono_cal_date_approx_cumulative_fct for approximation of the posterior cumulative distribution function (CDF) for a single calibrated radiocarbon date.

This function returns a vectorized CDF: it takes an array of scaled dates and returns an array of the same shape containing the corresponding cumulative probabilities.

The CDF is constructed from either:
- a callable density function density, or
- precomputed subdivision components (interval bounds, middle points, densities at middle points).

As in the scalar version, the continuity of the CDF is a direct mathematical consequence of the sampling strategy: intervals are chosen with probability proportional to their middle–point density, and values inside each interval are drawn uniformly. This implies a linear increase of the CDF inside each interval.

Parameters:

Name	Type	Description	Default
density	callable	A function density(dates: np.ndarray) -> np.ndarray computing the approximate posterior density at given points. Required if subdivision_components is None.	None
nb_intervals	int	Number of subintervals to discretize the support for the density approximation. Overrided internally if subdivision_components is provided.	1000
support_bounds	tuple of float	Tuple (min, max) specifying the bounds of the scaled date support. Only (0,1) is currently supported. Default is (0,1).	(0, 1)
subdivision_components	tuple or list of three numpy.ndarray	Precomputed components for the density approximation: - array of interval bounds (length N+1) - array of middle points (length N) - array of density values at middle points (length N) If provided, density is ignored and nb_intervals is overrided.	None

Returns:

Name	Type	Description
cumulative_density	callable	A vectorized function: CDF(dates: np.ndarray) -> np.ndarray returning the approximate posterior cumulative probabilities for the given dates.

Raises:

Type	Description
NotImplementedError	If support_bounds is not (0, 1).
ValueError	If neither density nor subdivision_components are provided, or if subdivision_components does not contain exactly three arrays.

Notes

Let the support \([a,b]\) be subdivided into N intervals \([x_j, x_{j+1}]\) with middle points \(m_j = \frac{x_j + x_{j+1}}{2}\), and let \(f_j = density(m_j)\). The approximate \(CDF\) at a point \(d \in [x_j, x_{j+1}]\) is: $$ CDF(d) = \sum_{i=1}^{j-1} \dfrac{f_i}{\sum_{k=1}^N f_k} + \dfrac{f_j}{\sum_{k=1}^N f_k} \dfrac{d - x_j}{x_{j+1} - x_j} \, , $$
where the first term sums contributions from previous intervals, and the second term accounts for the uniform distribution inside the current interval.

The CDF is continuous because the sampling strategy draws uniformly inside the selected interval.

The vectorized function applies this formula simultaneously to an array of points.

Examples:

>>> density_fn = mono_cal_date_approx_density(
...     mesure=0.954,
...     lab_error=0.002,
...     bnn_model=my_trained_bnn,
...     nb_curves=200
... )
>>> cdf_fn = mono_cal_date_approx_vect_cumulative_fct(density=density_fn, nb_intervals=1000)
>>> scaled_dates = np.array([0.25, 0.5, 0.8])
>>> cdf_values = cdf_fn(scaled_dates)
>>> cdf_values.shape
(3,)

Source code in src/bnn_for_14C_calibration/calibration_utils.py

def mono_cal_date_approx_vect_cumulative_fct(
    density: Optional[Callable[[np.ndarray], np.ndarray]] = None,
    nb_intervals: int = 1000,
    support_bounds: Tuple[float, float] = (0, 1),
    subdivision_components: Optional[Union[Tuple[np.ndarray, np.ndarray, np.ndarray], List[np.ndarray]]] = None
) -> Callable[[np.ndarray], np.ndarray]:
    """
    Vectorized version of `mono_cal_date_approx_cumulative_fct` for approximation 
    of the posterior cumulative distribution function (CDF)
    for a single calibrated radiocarbon date.

    This function returns a *vectorized* CDF: it takes an array of scaled dates
    and returns an array of the same shape containing the corresponding cumulative
    probabilities.

    The CDF is constructed from either:  
      - a callable density function `density`, or  
      - precomputed subdivision components (interval bounds, middle points,
        densities at middle points).

    As in the scalar version, the continuity of the CDF is a direct mathematical
    consequence of the sampling strategy: intervals are chosen with probability
    proportional to their middle–point density, and values inside each interval
    are drawn uniformly. This implies a linear increase of the CDF inside each
    interval.

    Parameters
    ----------
    density : callable, optional
        A function `density(dates: np.ndarray) -> np.ndarray` computing the 
        approximate posterior density at given points. Required if `subdivision_components` is None.
    nb_intervals : int, default=1000
        Number of subintervals to discretize the support for the density approximation.
        Overrided internally if `subdivision_components` is provided.
    support_bounds : tuple of float, optional
        Tuple `(min, max)` specifying the bounds of the scaled date support. Only `(0,1)` is
        currently supported. Default is `(0,1)`.
    subdivision_components : tuple or list of three numpy.ndarray, optional
        Precomputed components for the density approximation:  
            - array of interval bounds (length N+1)  
            - array of middle points  (length N)  
            - array of density values at middle points  (length N)  
        If provided, `density` is ignored and `nb_intervals` is overrided.

    Returns
    -------
    cumulative_density : callable
        A vectorized function:
            `CDF(dates: np.ndarray) -> np.ndarray` returning the approximate
            posterior cumulative probabilities for the given dates.

    Raises
    ------
    NotImplementedError
        If `support_bounds` is not `(0, 1)`.
    ValueError
        If neither `density` nor `subdivision_components` are provided, or if
        `subdivision_components` does not contain exactly three arrays.

    Notes
    -----
    Let the support $[a,b]$ be subdivided into N intervals $[x_j, x_{j+1}]$ with
    middle points $m_j = \\frac{x_j + x_{j+1}}{2}$, and let $f_j = density(m_j)$.
    The approximate $CDF$ at a point $d \in [x_j, x_{j+1}]$ is:
    $$
    CDF(d) = \sum_{i=1}^{j-1} \\dfrac{f_i}{\sum_{k=1}^N f_k}  + 
    \\dfrac{f_j}{\sum_{k=1}^N f_k} \\dfrac{d - x_j}{x_{j+1} - x_j} \\, ,
    $$  
    where the first term sums contributions from previous intervals, and the second
    term accounts for the uniform distribution inside the current interval. 

    The CDF is continuous because the sampling strategy draws uniformly inside
    the selected interval.

    The vectorized function applies this formula simultaneously to an array of
    points.  

    Examples
    --------
    >>> density_fn = mono_cal_date_approx_density(
    ...     mesure=0.954,
    ...     lab_error=0.002,
    ...     bnn_model=my_trained_bnn,
    ...     nb_curves=200
    ... )
    >>> cdf_fn = mono_cal_date_approx_vect_cumulative_fct(density=density_fn, nb_intervals=1000)
    >>> scaled_dates = np.array([0.25, 0.5, 0.8])
    >>> cdf_values = cdf_fn(scaled_dates)
    >>> cdf_values.shape
    (3,)

    """  

     # traitement du support et contrôle des arguments fournis
    if support_bounds != (0,1) :
        raise NotImplementedError(
            "Only the default support (0,1) is currently supported"
        )

    if density == None and subdivision_components == None :
        raise ValueError(
            "At least one of 'density' or 'subdivision_components' must be provided (not None)"
        )

    if subdivision_components != None and len(subdivision_components) != 3 :
        raise ValueError(
            "'subdivision_components' must be a tuple or list of length 3 containing arrays: interval bounds, middle points, and densities at middle points"
        )

    if subdivision_components != None : 

        intervals_bounds = subdivision_components[0]
        middle_points = subdivision_components[1]
        middle_points_density = subdivision_components[2]
        nb_intervals = len(middle_points)

    else :

        # subdivision du support en nb_intervals : calcul des bornes des sous-intervalles
        intervals_bounds = np.linspace(support_bounds[0], support_bounds[1], nb_intervals+1, dtype=np.float64)

        # évaluation de la densité aux points milieu
        # midle_points = (support_bounds[1] - support_bounds[0])/(2*nb_intervals) + intervals_bounds[:-1]
        middle_points = (intervals_bounds[:-1] + intervals_bounds[1:])/2
        middle_points_density = density(middle_points)

    # la borne inférieure de l'intervalle qui contient la date d sur laquelle évaluer la fonction
    lower_test = lambda d : intervals_bounds[:-1] <= d.reshape((-1,1)) 

    def idx(d) : 
        condition_res = np.where(lower_test(d))
        idx_of_idx = []
        for i in range(len(d)) :
            idx_of_idx.append(np.where(condition_res[0]==i)[0][-1]) #-1 car c'est la dernière borne pour laquelle le test ci-dessus vaut true
            # NB : np.where renvoie un tuple d'array (ici c'est un tuple avec un seul array) et c'est le premier élément du tuple qui nous intéresse ici
        idx_of_idx = np.array(idx_of_idx)

        res_idx = condition_res[1][idx_of_idx]
        return res_idx

    # # autre possibilité de calcul de la fonction idx sans recourir à une boucle :

    # # test sur les bornes inf et sup des intervalles qui contiennent les dates
    # lower_test = lambda d : intervals_bounds[:-1] <= d.reshape((-1,1))
    # intervals_bounds[-1] = intervals_bounds[-1] + 1 # on ajoute un nb positif à la borne sup du support pour donner un sens a inf <= d < sup pour tout d du support (sinon d = sup poserait problème)
    # upper_test = lambda d : d.reshape((-1,1)) < intervals_bounds[1:]

    # # enfin on trouverait les indices des bornes inf des intervalles comme suit :
    # idx = lambda d : np.where(lower_test(d)*upper_test(d))[1]
    # # NB : np.where renvoie un tuple d'array (ici c'est un tuple avec un 2 arrays : array des dimensions et array des indices) et c'est le deuxième élément du tuple qui nous intéresse ici

    # calcul de la densité cumulée au point d 
    h = (support_bounds[1] - support_bounds[0])/nb_intervals

    def cumulative_density(d):
        res_idx = idx(d)
        inf_cumulative_density = []
        for idx_d in res_idx :
            inf_cumulative_density.append(middle_points_density[:idx_d].sum())
        inf_cumulative_density = np.array(inf_cumulative_density)

        vect_cumulative_density = (middle_points_density[res_idx]*(d - intervals_bounds[res_idx])/h + inf_cumulative_density)/middle_points_density.sum()
        return vect_cumulative_density

    return cumulative_density

mono_cal_date_discrete_approx_quantile_fct(density=None, nb_intervals=1000, support_bounds=(0, 1), subdivision_components=None)

Approximate the posterior quantile function for a single calibrated radiocarbon date using a discrete grid derived from the posterior density approximation.

The quantile returned for a level alpha is the closest discrete realization available on the grid, i.e. the first grid point whose discretized CDF is greater than or equal to alpha. This is a purely discrete approximation, unlike the continuous quantile obtainable via inversion of a continuous CDF.

Parameters:

Name	Type	Description	Default
density	callable	A function density(dates: np.ndarray) -> np.ndarray computing the approximate posterior density at given points. Required if subdivision_components is None.	None
nb_intervals	int	Number of subintervals to discretize the support for the density approximation. Ignored if subdivision_components is provided.	1000
support_bounds	tuple of float	Tuple (min, max) specifying the bounds of the scaled date support. Only (0,1) is currently supported. Default is (0,1).	(0, 1)
subdivision_components	tuple or list of three numpy.ndarray	Precomputed components for the density approximation: - array of interval bounds (length N+1) - array of middle points (length N) - array of density values at middle points (length N) If provided, density is ignored.	None

Returns:

Name	Type	Description
discrete_alpha_quantile	callable	A function Q(alpha: float) -> float returning the discrete approximation of the posterior quantile of order alpha.

Raises:

Type	Description
NotImplementedError	If support_bounds is not (0, 1).
ValueError	If neither density nor subdivision_components are provided, or if subdivision_components does not contain exactly three arrays.

Notes

Approximation of the CDF.
Let the posterior density be approximated on nb_intervals subintervals by evaluating it at middle points \(m_j\) with corresponding values \(f_j\). The cumulative density assigned to a middle point is then written as: $$ \forall j \in \{ 1, \cdots, N \} \,, F(m_j) = \dfrac{\frac{f_j}{2} + \sum_{i=1}^{j-1} f_i}{\sum_{k=1}^N f_k}, $$

where:
- \(\sum_{i<j} f_i\) represents the (unnormalized) cumulative probability mass of the intervals preceding the j-th,
- \(f_j/2\) accounts for integrating half of the (piecewise) constant density over the current interval up to its midpoint.

This expression arises directly from the integration of the piecewise constant density implicitly defined by the sampling strategy of the approximate posterior: selecting an interval proportionally to its density and drawing uniformly inside it.

Discrete quantile.
Extending the discretized CDF with values:
- \(F = 0\) at the lower bound,
- \(F = 1\) at the upper bound,
we obtain an ordered set of cumulative probabilities: $$ 0 = F_0 < F(m_1) < \dots < F(m_{N}) < F_{N+1} = 1. $$

The discrete quantile is then defined as: $$ Q_{\text{disc}}(\alpha) = x_{k} \quad \text{where } k = \min\{ i : F_i \ge \alpha \}, $$ with \(x_k\) the corresponding support point (a midpoint, the lower bound, or the upper bound).

This approximation converges to the true quantile as the number of intervals increases, because the piecewise-constant density and its cumulative sum are standard Riemann approximations of the continuous density and CDF.

Examples:

>>> density_fn = mono_cal_date_approx_density(
...     mesure=0.954,
...     lab_error=0.002,
...     bnn_model=my_trained_bnn,
...     nb_curves=200
... )
>>> q_fn = mono_cal_date_discrete_approx_quantile_fct(density_fn, nb_intervals=2000)
>>> q25 = q_fn(0.25)
>>> isinstance(q25, float)
True

Source code in src/bnn_for_14C_calibration/calibration_utils.py

def mono_cal_date_discrete_approx_quantile_fct(
    density: Optional[Callable[[np.ndarray], np.ndarray]] = None,
    nb_intervals: int = 1000,
    support_bounds: Tuple[float, float] = (0, 1),
    subdivision_components: Optional[
        Union[
            Tuple[np.ndarray, np.ndarray, np.ndarray],
            List[np.ndarray]
        ]
    ] = None
) -> Callable[[float], float]:
    """
    Approximate the posterior quantile function for a single calibrated radiocarbon date
    using a discrete grid derived from the posterior density approximation.

    The quantile returned for a level ``alpha`` is the *closest discrete realization*
    available on the grid, i.e. the first grid point whose discretized CDF is greater
    than or equal to ``alpha``. This is a purely discrete approximation, unlike the 
    continuous quantile obtainable via inversion of a continuous CDF.

    Parameters
    ----------
    density : callable, optional
        A function `density(dates: np.ndarray) -> np.ndarray` computing the 
        approximate posterior density at given points. Required if `subdivision_components` is None.
    nb_intervals : int, default=1000
        Number of subintervals to discretize the support for the density approximation.
        Ignored if `subdivision_components` is provided.
    support_bounds : tuple of float, optional
        Tuple `(min, max)` specifying the bounds of the scaled date support. Only `(0,1)` is
        currently supported. Default is `(0,1)`.
    subdivision_components : tuple or list of three numpy.ndarray, optional
        Precomputed components for the density approximation:  
            - array of interval bounds (length N+1)  
            - array of middle points  (length N)  
            - array of density values at middle points  (length N)  
        If provided, `density` is ignored.

    Returns
    -------
    discrete_alpha_quantile : callable
        A function `Q(alpha: float) -> float` returning the discrete approximation
        of the posterior quantile of order ``alpha``.

    Raises
    ------
    NotImplementedError
        If `support_bounds` is not `(0, 1)`.
    ValueError
        If neither `density` nor `subdivision_components` are provided, or if
        `subdivision_components` does not contain exactly three arrays.

    Notes
    -----
    **Approximation of the CDF.**  
    Let the posterior density be approximated on ``nb_intervals`` subintervals
    by evaluating it at middle points $m_j$ with corresponding values $f_j$. 
    The cumulative density assigned to a middle point is then written as:
    $$
        \\forall j \in \\\{ 1, \cdots, N \\\} \,,
         F(m_j)
        = \\dfrac{\\frac{f_j}{2} + \sum_{i=1}^{j-1} f_i}{\sum_{k=1}^N f_k},
    $$

    where:  
    - $\sum_{i<j} f_i$ represents the (unnormalized) cumulative probability mass of the intervals 
      preceding the j-th,  
    - $f_j/2$ accounts for integrating half of the (piecewise) constant density over the current 
      interval up to its midpoint.

    This expression **arises directly** from the integration of the *piecewise constant*
    density implicitly defined by the sampling strategy of the approximate posterior:
    selecting an interval proportionally to its density and drawing uniformly inside it.

    **Discrete quantile.**  
    Extending the discretized CDF with values:  
        - $F = 0$ at the lower bound,  
        - $F = 1$ at the upper bound,  
    we obtain an ordered set of cumulative probabilities:
    $$
    0 = F_0 < F(m_1) < \dots < F(m_{N}) < F_{N+1} = 1.
    $$

    The discrete quantile is then defined as:
    $$
        Q_{\\text{disc}}(\\alpha)
        = x_{k} \quad \\text{where } k = \min\\\{ i : F_i \ge \\alpha \\\},
    $$
    with $x_k$ the corresponding support point (a midpoint, the lower bound, 
    or the upper bound).

    This approximation converges to the true quantile as the number of intervals 
    increases, because the piecewise-constant density and its cumulative sum 
    are standard Riemann approximations of the continuous density and CDF.

    Examples
    --------
    >>> density_fn = mono_cal_date_approx_density(
    ...     mesure=0.954,
    ...     lab_error=0.002,
    ...     bnn_model=my_trained_bnn,
    ...     nb_curves=200
    ... )
    >>> q_fn = mono_cal_date_discrete_approx_quantile_fct(density_fn, nb_intervals=2000)
    >>> q25 = q_fn(0.25)
    >>> isinstance(q25, float)
    True
    """

    # traitement du support et contrôle des arguments fournis
    if support_bounds != (0,1) :
        raise NotImplementedError(
            "Only the default support (0,1) is currently supported"
        )

    if density == None and subdivision_components == None :
        raise ValueError(
            "At least one of 'density' or 'subdivision_components' must be provided (not None)"
        )

    if subdivision_components != None and len(subdivision_components) != 3 :
        raise ValueError(
            "'subdivision_components' must be a tuple or list of length 3 containing arrays: interval bounds, middle points, and densities at middle points"
        )

    if subdivision_components != None : 

        intervals_bounds = subdivision_components[0]
        middle_points = subdivision_components[1]
        middle_points_density = subdivision_components[2]
        # nb_intervals = len(middle_points)

    else :

        # subdivision du support en nb_intervals : calcul des bornes des sous-intervalles
        intervals_bounds = np.linspace(support_bounds[0], support_bounds[1], nb_intervals+1, dtype=np.float64)

        # évaluation de la densité aux points milieu
        # midle_points = (support_bounds[1] - support_bounds[0])/(2*nb_intervals) + intervals_bounds[:-1]
        middle_points = (intervals_bounds[:-1] + intervals_bounds[1:])/2
        middle_points_density = density(middle_points)

    # calcul de la fonction de répartition aux points milieu
    middle_points_density_cumsum = middle_points_density.cumsum()
    middle_points_cumulative_density = (middle_points_density/2 + np.concatenate((np.array([0]), middle_points_density_cumsum[:-1]), dtype=np.float64)) / middle_points_density.sum()

    # Enfin on peut déterminer les fonctions de répartition et quantile ainsi discrétisées
    discrete_cumulative_density = np.concatenate((np.array([0]), middle_points_cumulative_density, np.array([1])), dtype=np.float64)
    discrete_quantiles = np.concatenate((np.array([support_bounds[0]]), middle_points, np.array([support_bounds[1]])), dtype=np.float64)

    # Pour terminer, on définit la fonction quantile qui retournera le quantile d'ordre alpha sur base de cette discrétisation
    idx = lambda alpha : np.where(discrete_cumulative_density >= alpha)[0][0] # les densités cumulées et les quantiles étant ordonnés, il suffit de prendre le premier indice (autrement il aurait fallu ordonner d'abord)
    discrete_alpha_quantile = lambda alpha : discrete_quantiles[idx(alpha)]

    return discrete_alpha_quantile

mono_cal_date_exact_approx_quantile_fct(density=None, nb_intervals=1000, support_bounds=(0, 1), subdivision_components=None)

Compute a continuous approximation of the posterior quantile function for a single calibrated radiocarbon date by analytically inverting the piecewise-linear CDF obtained from the piecewise-constant approximation of the posterior density.

This contrasts with the discrete quantile approximation, which selects the closest grid point. Here, the quantile is obtained by solving an affine equation inside the unique interval where the continuous CDF crosses alpha.

Parameters:

Name	Type	Description	Default
density	callable	A function density(dates: np.ndarray) -> np.ndarray computing the approximate posterior density at given points. Required if subdivision_components is None.	None
nb_intervals	int	Number of subintervals to discretize the support for the density approximation. Overrided internally if subdivision_components is provided.	1000
support_bounds	tuple of float	Tuple (min, max) specifying the bounds of the scaled date support. Only (0,1) is currently supported. Default is (0,1).	(0, 1)
subdivision_components	tuple or list of three numpy.ndarray	Precomputed components for the density approximation: - array of interval bounds (length N+1) - array of middle points (length N) - array of density values at middle points (length N) If provided, density is ignored and nb_intervals is overrided.	None

Returns:

Name	Type	Description
exact_alpha_quantile	callable	A function Q(alpha: float) -> float returning the continuous approximate posterior quantile of order alpha.

Raises:

Type	Description
NotImplementedError	If support_bounds is not (0, 1).
ValueError	If neither density nor subdivision_components are provided, or if subdivision_components does not contain exactly three arrays.

Notes

1. Approximate posterior density.
As in the discrete quantile version, the posterior density is approximated by a piecewise-constant function: $$ f(d) \approx f(m_j) := f_j \quad \text{for } d \in [x_j, x_{j+1}], $$ where \(x_j\) are interval bounds and \(m_j = \frac{x_j+x_{j+1}}{2}\) are midpoints.

2. Approximate CDF at interval bounds.
The cumulative probability at the bounds is: $$ F(x_j) = \dfrac{\sum_{i=1}^{j} f_i}{\sum_{k=1}^{N} f_k}. $$ These values form a strictly increasing sequence from 0 to 1.

3. Continuous quantile: inversion on each interval.
Inside the interval \([x_j, x_{j+1}]\), the density is constant, hence the CDF is affine: $$ F(d) = F(x_j) + \dfrac{f_j}{\sum_{k=1}^N f_k}\,\dfrac{d - x_j}{x_{j+1} - x_j}. $$ To find the quantile of order \(\alpha\), determine the unique interval where: $$ F(x_j) < \alpha \le F(x_{j+1}), $$ and solve for d: $$ d = x_j + \frac{h}{f_j} \left( \alpha \sum_{k=1}^N f_k- \sum_{i=1}^{j} f_i \right) := Q(\alpha) \, , $$ where \(h = x_{j+1} - x_j \, , \quad \forall j \in \{ 1, \cdots, N \}\).

4. Difference with the discrete quantile.
- Discrete quantile selects the closest grid point.
- Continuous quantile solves a linear equation inside the interval.
- This yields a true continuous function of \(\alpha\) with no jumps.
- Both converge to the exact quantile as the grid is refined, but the continuous version converges faster.

Examples:

>>> density_fn = mono_cal_date_approx_density(
...     mesure=0.954,
...     lab_error=0.002,
...     bnn_model=my_trained_bnn,
...     nb_curves=200
... )
>>> q_fn = mono_cal_date_exact_approx_quantile_fct(density_fn, nb_intervals=2000)
>>> q25 = q_fn(0.25)
>>> isinstance(q25, float)
True

Source code in src/bnn_for_14C_calibration/calibration_utils.py

def mono_cal_date_exact_approx_quantile_fct(
    density: Optional[Callable[[np.ndarray], np.ndarray]] = None,
    nb_intervals: int = 1000,
    support_bounds: Tuple[float, float] = (0, 1),
    subdivision_components: Optional[
        Union[
            Tuple[np.ndarray, np.ndarray, np.ndarray],
            List[np.ndarray]
        ]
    ] = None
) -> Callable[[float], float]:
    """
    Compute a *continuous* approximation of the posterior quantile function 
    for a single calibrated radiocarbon date by analytically inverting the
    piecewise-linear CDF obtained from the piecewise-constant approximation 
    of the posterior density.

    This contrasts with the *discrete* quantile approximation, which selects
    the closest grid point. Here, the quantile is obtained by *solving an 
    affine equation inside the unique interval where the continuous CDF 
    crosses* ``alpha``.

    Parameters
    ----------
    density : callable, optional
        A function `density(dates: np.ndarray) -> np.ndarray` computing the 
        approximate posterior density at given points. Required if `subdivision_components` is None.
    nb_intervals : int, default=1000
        Number of subintervals to discretize the support for the density approximation.
        Overrided internally if `subdivision_components` is provided.
    support_bounds : tuple of float, optional
        Tuple `(min, max)` specifying the bounds of the scaled date support. Only `(0,1)` is
        currently supported. Default is `(0,1)`.
    subdivision_components : tuple or list of three numpy.ndarray, optional
        Precomputed components for the density approximation:  
            - array of interval bounds (length N+1)  
            - array of middle points  (length N)  
            - array of density values at middle points  (length N)  
        If provided, `density` is ignored and `nb_intervals` is overrided.

    Returns
    -------
    exact_alpha_quantile : callable  
        A function `Q(alpha: float) -> float` returning the *continuous* 
        approximate posterior quantile of order ``alpha``.

    Raises
    ------
    NotImplementedError
        If `support_bounds` is not `(0, 1)`.
    ValueError
        If neither `density` nor `subdivision_components` are provided, or if
        `subdivision_components` does not contain exactly three arrays.

    Notes
    -----
    **1. Approximate posterior density.**  
    As in the discrete quantile version, the posterior density is approximated
    by a piecewise-constant function:
    $$
        f(d) \\approx f(m_j) := f_j \quad \\text{for } d \in [x_j, x_{j+1}],
    $$
    where $x_j$ are interval bounds and $m_j = \\frac{x_j+x_{j+1}}{2}$ are midpoints.

    **2. Approximate CDF at interval bounds.**  
    The cumulative probability at the bounds is:
    $$
        F(x_j)
        = \\dfrac{\sum_{i=1}^{j} f_i}{\sum_{k=1}^{N} f_k}.
    $$
    These values form a strictly increasing sequence from 0 to 1.

    **3. Continuous quantile: inversion on each interval.**  
    Inside the interval $[x_j, x_{j+1}]$, the density is constant, hence
    the CDF is affine:
    $$
        F(d) = F(x_j) + \\dfrac{f_j}{\sum_{k=1}^N f_k}\,\\dfrac{d - x_j}{x_{j+1} - x_j}.
    $$
    To find the quantile of order $\\alpha$, determine the unique interval where:
    $$
        F(x_j) < \\alpha \le F(x_{j+1}),
    $$
    and solve for d:
    $$
        d = x_j + 
            \\frac{h}{f_j} \left( \\alpha \sum_{k=1}^N f_k- \sum_{i=1}^{j} f_i \\right)
            := Q(\\alpha) \, ,
    $$
    where $h = x_{j+1} - x_j \, , \quad \\forall j \in \\{ 1, \cdots, N \\}$.

    **4. Difference with the discrete quantile.**  
    - *Discrete quantile* selects the **closest grid point**.  
    - *Continuous quantile* solves a **linear equation** inside the interval.  
    - This yields a *true continuous function* of $\\alpha$ with no jumps.  
    - Both converge to the exact quantile as the grid is refined, but the
      continuous version converges faster.

    Examples
    --------
    >>> density_fn = mono_cal_date_approx_density(
    ...     mesure=0.954,
    ...     lab_error=0.002,
    ...     bnn_model=my_trained_bnn,
    ...     nb_curves=200
    ... )
    >>> q_fn = mono_cal_date_exact_approx_quantile_fct(density_fn, nb_intervals=2000)
    >>> q25 = q_fn(0.25)
    >>> isinstance(q25, float)
    True
    """

    # traitement du support et contrôle des arguments fournis
    if support_bounds != (0,1) :
        raise NotImplementedError(
            "Only the default support (0,1) is currently supported"
        )

    if density == None and subdivision_components == None :
        raise ValueError(
            "At least one of 'density' or 'subdivision_components' must be provided (not None)"
        )

    if subdivision_components != None and len(subdivision_components) != 3 :
        raise ValueError(
            "'subdivision_components' must be a tuple or list of length 3 containing arrays: interval bounds, middle points, and densities at middle points"
        )

    if subdivision_components != None : 

        intervals_bounds = subdivision_components[0]
        middle_points = subdivision_components[1]
        middle_points_density = subdivision_components[2]
        nb_intervals = len(middle_points)

    else :

        # subdivision du support en nb_intervals : calcul des bornes des sous-intervalles
        intervals_bounds = np.linspace(support_bounds[0], support_bounds[1], nb_intervals+1, dtype=np.float64)

        # évaluation de la densité aux points milieu
        # midle_points = (support_bounds[1] - support_bounds[0])/(2*nb_intervals) + intervals_bounds[:-1]
        middle_points = (intervals_bounds[:-1] + intervals_bounds[1:])/2
        middle_points_density = density(middle_points)

    # calcul de la fonction de répartition aux "nb_intervals + 1" bornes des sous-intervalles du support
    middle_points_density_sum = middle_points_density.sum()
    middle_points_density_cumsum = middle_points_density.cumsum()
    intervals_bounds_cumulative_density = np.concatenate((np.array([0]), middle_points_density_cumsum / middle_points_density_sum), dtype=np.float64)

    # on met à 1 les éventuelles valeurs de la densité cumulée (fonction de répartition) qui seraient supérieures à 1 à cause des erreurs d'arrondis
    intervals_bounds_cumulative_density = np.where(intervals_bounds_cumulative_density > 1., 1., intervals_bounds_cumulative_density)

    # on rajoute 0 comme premier élément de 'middle_points_density_cumsum' pour que cela puisse bien marcher avec les indices lors du calcul de 'd_alpha' plus loin
    # en effet on voudra que si idx_borne_inf = 0 (donc la borne inf de l'intervalle = borne inf du support), alors middle_points_density_cumsum[idx_borne_inf] doit donner 0
    middle_points_density_cumsum = np.concatenate((np.array([0]), middle_points_density_cumsum), dtype=np.float64)
    # (attention : len(middle_points_density_cumsum) est nb_intervals + 1 désormais alors que len(middle_points_density) vaut toujours nb_intervals, ce qui fait en sorte que 
    # middle_points_density[idx_borne_inf] donne bien la densité (non normalisée) au point milieu de l'intervalle de borne inf = intervals_bounds[idx_borne_inf], ce qui est 
    # exactement le résultat voulus)

    # on peut enfin calculer notre fonction quantile
    def exact_alpha_quantile(alpha) :
        try :
            # Si l'une des "nb_intervals + 1" bornes a une densité cumulée (fonction de répartition) qui vaut alpha, alors c'est le quantile recherché
            idx_borne = np.where(intervals_bounds_cumulative_density == alpha)[0][0]
            return intervals_bounds[idx_borne]
        except IndexError :
            # Sinon, aucune des bornes n'est le quantile :
            # on cherche alors la borne inf de l'intervalle qui contient notre quantile : c'est la dernère borne avec une densité cumulée < alpha
            idx_borne_inf = np.where(intervals_bounds_cumulative_density < alpha)[0][-1]

            # on calcule alors le quantile d'ordre alpha, d_alpha, suivant la formule de la fonction quantile 
            h = (support_bounds[1] - support_bounds[0])/(nb_intervals)
            d_alpha = intervals_bounds[idx_borne_inf] + h / middle_points_density[idx_borne_inf] * (alpha * middle_points_density_sum - middle_points_density_cumsum[idx_borne_inf])
            return d_alpha

    # Enfin on peut retourner le fonction quantile ainsi construite
    return exact_alpha_quantile

multi_cal_date_approx_density(mesures, lab_errors, bnn_model, nb_curves=100, prior_density='default', batch_size=None)

Approximate the joint posterior density for multiple radiocarbon dates.

This function generalizes the single-date calibration approach to the case where several radiocarbon dates are calibrated simultaneously.
It computes an approximate joint posterior density over the vector of (scaled) calendar dates, using a trained Bayesian Neural Network (BNN) model as an estimator of the calibration curve in the F¹⁴C domain.

Parameters:

Name	Type	Description	Default
mesures	np.ndarray of shape (n_dates,)	The measured radiocarbon ages expressed in the F¹⁴C domain.	required
lab_errors	np.ndarray of shape (n_dates,)	The laboratory measurement uncertainties (standard deviations), also expressed in the F¹⁴C domain.	required
bnn_model	object	The trained Bayesian Neural Network model used to estimate the predictive distribution. Must be compatible with bnn_make_predictions_ and approximate the calibration curve in the F¹⁴C domain.	required
nb_curves	int	Number of stochastic realizations (Monte Carlo samples) to use for approximating the BNN predictive distribution. Default is 100.	100
prior_density	(default, callable)	Prior probability density over the vector of scaled calendar dates. - "default": a uniform prior over the hypercube [0, 1]^n_dates. - callable: a custom prior density function of the form f(dates) → np.ndarray. Default is "default".	"default"
batch_size	int	Batch size for model predictions, passed to the internal bnn_make_predictions_ function. Default is None.	None

Returns:

Name	Type	Description
density	callable	A function density(dates: np.ndarray) -> np.ndarray that computes the (unnormalized) approximate joint posterior density of calibrated dates. Each row in dates corresponds to one vector of scaled calendar dates. The density is known up to a normalization constant.

Notes

• The posterior density is proportional to:
  \( p(\mathbf{d}|\mathbf{m}) ∝ p(\mathbf{d}) × E_{BNN}[ \prod_i \exp(-(m_i - \hat{F}^{14}C(d_i))^2 / (2σ_i^2)) ] \).
  The expectation over the BNN distribution is approximated by Monte Carlo averaging.
• The "default" prior corresponds to a uniform independent prior over each scaled date in [0, 1].
• The output density is not normalized; handling normalization here via numerical integration over a multi-dimensional grid is not a good idea because of the curse of dimensionality. In general, trying to do so will not be necessary because the outputs densities are expected to be used within a MCMC sampler (e.g. Metropolis-Hastings within Gibbs sampler).

Raises:

Type	Description
NotImplementedError	If a non-default prior density is provided (custom priors are not yet supported).

Examples:

>>> mesures = np.array([0.954, 0.928])
>>> lab_errors = np.array([0.002, 0.003])
>>> density_fn = multi_cal_date_approx_density(
...     mesures=mesures,
...     lab_errors=lab_errors,
...     bnn_model=my_trained_bnn,
...     nb_curves=200
... )
>>> date_grid = np.random.rand(100, 2)  # 100 candidate date vectors
>>> posterior_vals = density_fn(date_grid)
>>> posterior_vals.shape
(100,)

Source code in src/bnn_for_14C_calibration/calibration_utils.py

def multi_cal_date_approx_density(
    mesures: np.ndarray,
    lab_errors: np.ndarray,
    bnn_model: object,
    nb_curves: int = 100,
    prior_density: Union[str, Callable[[np.ndarray], np.ndarray]] = "default",
    batch_size: Optional[int] = None
) -> Callable[[np.ndarray], np.ndarray]:
    """
    Approximate the joint posterior density for multiple radiocarbon dates.

    This function generalizes the single-date calibration approach to the case 
    where several radiocarbon dates are calibrated simultaneously.  
    It computes an approximate joint posterior density over the vector of 
    (scaled) calendar dates, using a trained Bayesian Neural Network (BNN) model 
    as an estimator of the calibration curve in the F¹⁴C domain.

    Parameters
    ----------
    mesures : np.ndarray of shape (n_dates,)
        The measured radiocarbon ages expressed in the F¹⁴C domain.
    lab_errors : np.ndarray of shape (n_dates,)
        The laboratory measurement uncertainties (standard deviations), 
        also expressed in the F¹⁴C domain.
    bnn_model : object
        The trained Bayesian Neural Network model used to estimate the predictive distribution.  
        Must be compatible with `bnn_make_predictions_` and approximate the calibration curve in the F¹⁴C domain.
    nb_curves : int, optional
        Number of stochastic realizations (Monte Carlo samples) to use for 
        approximating the BNN predictive distribution.  
        Default is `100`.
    prior_density : {"default", callable}, optional
        Prior probability density over the vector of scaled calendar dates.  
        - `"default"`: a uniform prior over the hypercube `[0, 1]^n_dates`.  
        - `callable`: a custom prior density function of the form `f(dates) → np.ndarray`.  
        Default is `"default"`.
    batch_size : int, optional
        Batch size for model predictions, passed to the internal 
        `bnn_make_predictions_` function.  
        Default is `None`.

    Returns
    -------
    density : callable
        A function `density(dates: np.ndarray) -> np.ndarray` that computes 
        the (unnormalized) approximate **joint posterior density** of calibrated dates.  
        Each row in `dates` corresponds to one vector of scaled calendar dates.  
        The density is known up to a normalization constant.

    Notes
    -----
    - The posterior density is proportional to:  
      \\( p(\\mathbf{d}|\\mathbf{m}) ∝ p(\\mathbf{d}) × E_{BNN}[ \\prod_i \\exp(-(m_i - \\hat{F}^{14}C(d_i))^2 / (2σ_i^2)) ] \\).  
      The expectation over the BNN distribution is approximated by Monte Carlo averaging.
    - The `"default"` prior corresponds to a uniform independent prior over each scaled date in `[0, 1]`.
    - The output density is **not normalized**; handling normalization here via numerical integration over a 
        multi-dimensional grid is not a good idea because of the curse of dimensionality. In general, trying 
        to do so will not be necessary because the outputs densities are expected to be used within 
        a MCMC sampler (e.g. Metropolis-Hastings within Gibbs sampler).

    Raises
    ------
    NotImplementedError
        If a non-default prior density is provided (custom priors are not yet supported).

    Examples
    --------
    >>> mesures = np.array([0.954, 0.928])
    >>> lab_errors = np.array([0.002, 0.003])
    >>> density_fn = multi_cal_date_approx_density(
    ...     mesures=mesures,
    ...     lab_errors=lab_errors,
    ...     bnn_model=my_trained_bnn,
    ...     nb_curves=200
    ... )
    >>> date_grid = np.random.rand(100, 2)  # 100 candidate date vectors
    >>> posterior_vals = density_fn(date_grid)
    >>> posterior_vals.shape
    (100,)
    """

    dim_dates = mesures.shape[0] # = len(mesures)
    # traitement de la densité à piori :
    if prior_density == "default" :
        support_lower_bound = np.array([0.] * dim_dates) # à remplacer par min_Xtrain ou min_Xtrain_val ou min_Xtrain_val_test plus tard suivant le cas ou date minimale gobale possible pour la calibration
        support_upper_bound = np.array([1.] * dim_dates) # à remplacer par max_Xtrain ou max_Xtrain_val ou max_Xtrain_val_test plus tard suivant le cas ou date maximale gobale possible pour la calibration
        prior_density = lambda d : np.float64((support_lower_bound <= d) * (d <= support_upper_bound)).prod(axis=1)/(support_upper_bound - support_lower_bound).prod()
    else :
        raise NotImplementedError(
            "Custom prior densities for multiple dates are not yet supported."
        )

    # predictions avec le modèle
    # d sera une matrice (un array 2-D numpy) dont chaque ligne correspond à un vecteur de dates sur lequel sera évaluée la densité jointe
    # predicted renvoie un array 2-D de taille (d.shape[0], dim_dates, nb_curves)
    predicted = lambda d : bnn_make_predictions_(bnn_model = bnn_model, X_test = d.reshape((-1,1)), iterations = nb_curves, batch_size = batch_size).reshape((-1,dim_dates,nb_curves))

    # densité approchée (connue à une constante près)
    mesures_broadcasted = mesures.repeat(nb_curves).reshape((-1,nb_curves))
    lab_errors_broadcasted = lab_errors.repeat(nb_curves).reshape(-1,nb_curves)

    density = lambda d : prior_density(d) * np.exp(-(mesures_broadcasted - predicted(d))**2/(2*lab_errors_broadcasted**2)).prod(axis=1, dtype=np.float64).mean(axis=1, dtype=np.float64) / (lab_errors.prod() * np.sqrt(2*np.pi)**dim_dates)

    return density

optimise_credible_interval(quantile, alpha)

Optimize the (1 - alpha)-credible interval for a single calibrated radiocarbon date by minimizing its length over all intervals of the form [Q(beta), Q(1 - alpha + beta)], where Q is a continuous posterior quantile function.

The optimization variable beta ∈ [0, alpha] determines the lower tail mass excluded from the interval. For a symmetric posterior, the optimal value is beta = alpha/2 (equal-tailed interval). For asymmetric posteriors, the optimal beta shifts to achieve the shortest interval of posterior mass (1 - alpha).

Parameters:

Name	Type	Description	Default
quantile	callable	A continuous quantile function Q(u) mapping u ∈ [0,1] to posterior dates. Must satisfy: - Q is non-decreasing, - Q(0) is the lower bound of the support, - Q(1) is the upper bound of the support.	required
alpha	float	Posterior tail probability. The credible interval contains mass (1 - alpha). Must satisfy 0 <= alpha <= 1.	required

Returns:

Name	Type	Description
beta_opt	object	The SciPy optimization result. The optimal value is: beta_opt.x[0] The corresponding credible interval is: [ Q(beta_opt.x[0]), Q(1 - alpha + beta_opt.x[0]) ].

Raises:

Type	Description
ValueError	If alpha is not between 0 and 1. If the quantile function does not accept valid inputs.

Notes

Mathematical background.

For any beta ∈ [0, alpha], the interval I_beta = [Q(beta), Q(1 - alpha + beta)] has posterior probability (1 - alpha).

Its length is L(beta) = Q(1 - alpha + beta) - Q(beta).

The shortest credible interval is obtained by solving: beta* = argmin_{beta ∈ [0, alpha]} L(beta).

When the posterior is unimodal, this interval coincides with the HPD (Highest Posterior Density) region. When the posterior is multimodal, the HPD region may be disconnected, but this function returns the shortest connected interval of mass (1 - alpha).

The optimization is performed using the Nelder–Mead method. In recent versions of SciPy, the bounds provided are used to constrain the shape and updates of the simplex, ensuring that the iterates remain within the allowed domain [0, alpha].

Examples:

>>> Q = lambda u: u**2      # toy quantile function
>>> res = optimise_credible_interval(Q, alpha=0.2)
>>> float(res.x[0]) >= 0
True

Source code in src/bnn_for_14C_calibration/calibration_utils.py

def optimise_credible_interval(
    quantile: Callable[[float], float],
    alpha: float
) -> object:
    """
    Optimize the (1 - alpha)-credible interval for a single calibrated radiocarbon
    date by minimizing its length over all intervals of the form
    [Q(beta), Q(1 - alpha + beta)], where Q is a continuous posterior quantile
    function.

    The optimization variable beta ∈ [0, alpha] determines the lower tail mass
    excluded from the interval. For a symmetric posterior, the optimal value is
    beta = alpha/2 (equal-tailed interval). For asymmetric posteriors, the optimal
    beta shifts to achieve the *shortest* interval of posterior mass (1 - alpha).

    Parameters
    ----------
    quantile : callable
        A continuous quantile function Q(u) mapping u ∈ [0,1] to posterior dates.
        Must satisfy:  
            - Q is non-decreasing,  
            - Q(0) is the lower bound of the support,  
            - Q(1) is the upper bound of the support.  
    alpha : float
        Posterior tail probability. The credible interval contains mass (1 - alpha).
        Must satisfy 0 <= alpha <= 1.

    Returns
    -------
    beta_opt : object
        The SciPy optimization result. The optimal value is:
            beta_opt.x[0]  
        The corresponding credible interval is:
            [ Q(beta_opt.x[0]), Q(1 - alpha + beta_opt.x[0]) ].

    Raises
    ------
    ValueError
        If alpha is not between 0 and 1.  
        If the quantile function does not accept valid inputs.

    Notes
    -----
    **Mathematical background.**

    For any beta ∈ [0, alpha], the interval
        I_beta = [Q(beta), Q(1 - alpha + beta)]
    has posterior probability (1 - alpha).

    Its length is
        L(beta) = Q(1 - alpha + beta) - Q(beta).

    The shortest credible interval is obtained by solving:
        beta* = argmin_{beta ∈ [0, alpha]} L(beta).

    When the posterior is unimodal, this interval coincides with the HPD
    (Highest Posterior Density) region. When the posterior is multimodal,
    the HPD region may be disconnected, but this function returns the shortest
    *connected* interval of mass (1 - alpha).

    The optimization is performed using the Nelder–Mead method. In recent
    versions of SciPy, the bounds provided are used to constrain the shape
    and updates of the simplex, ensuring that the iterates remain within the
    allowed domain [0, alpha].

    Examples
    --------
    >>> Q = lambda u: u**2      # toy quantile function
    >>> res = optimise_credible_interval(Q, alpha=0.2)
    >>> float(res.x[0]) >= 0
    True
    """

    # contrôle de alpha
    if not (0 <= alpha <= 1):
        raise ValueError("'alpha' must satisfy 0 <= alpha <= 1.")

    # longueur de l'intervalle de crédibilité pour un beta donné
    interval_length = lambda beta: quantile(1 - alpha + beta) - quantile(beta)

    # optimisation par Nelder-Mead (bornes gérées par la définition du simplexe)
    beta_opt = minimize(
        fun=interval_length,
        x0=np.array([alpha / 2]),
        method='Nelder-Mead',
        bounds=[(0., alpha)]
    )

    return beta_opt