diff --git a/src/pyuncertainnumber/calibration/tmcmc.py b/src/pyuncertainnumber/calibration/tmcmc.py
index ed0f9ae6..afaa0c97 100644
--- a/src/pyuncertainnumber/calibration/tmcmc.py
+++ b/src/pyuncertainnumber/calibration/tmcmc.py
@@ -1203,6 +1203,23 @@ def get_hdi_bounds(
     return pd.DataFrame(out).rename_axis("column").reset_index()
 
 
+def filter_hdi_bounds(df, level: int = 95) -> pd.DataFrame:
+    """Extract HDI bounds for a given credibility level.
+
+    args:
+        df (pd.DataFrame): HDI dataframe containing columns like `hdi_95_low`, `hdi_95_high`
+
+        level (int | str): significance level (e.g. 95, 90, 20)
+
+    returns:
+        (pd.DataFrame): DataFrame with columns [hdi_<level>_low, hdi_<level>_high]
+
+    """
+    level = str(level)
+    cols = [f"hdi_{level}_low", f"hdi_{level}_high"]
+    return df.loc[:, cols]
+
+
 def transform_old_trace_to_new(trace: list[list]) -> list[Stage]:
     """Transform old trace format to new Stage class format.
 
diff --git a/src/pyuncertainnumber/pba/aggregation.py b/src/pyuncertainnumber/pba/aggregation.py
index 7a49edaa..b7c5655c 100644
--- a/src/pyuncertainnumber/pba/aggregation.py
+++ b/src/pyuncertainnumber/pba/aggregation.py
@@ -135,6 +135,7 @@ def stochastic_mixture(
         return mixture_pbox(*converted_constructs, weights)
 
 
+# TODO. weights cannot be None for intervals with same weights
 def stacking(
     vec_interval: Interval | list[Interval],
     *,
diff --git a/src/pyuncertainnumber/pba/core.py b/src/pyuncertainnumber/pba/core.py
index e85c93cc..f2505a14 100644
--- a/src/pyuncertainnumber/pba/core.py
+++ b/src/pyuncertainnumber/pba/core.py
@@ -2,12 +2,14 @@
 from typing import TYPE_CHECKING
 
 from abc import ABC, abstractmethod
-import numpy as np
 from numpy.typing import ArrayLike
 import scipy.stats as sps
 from numbers import Number
 from pyuncertainnumber.pba.pbox_abc import Pbox, Staircase
+from pyuncertainnumber.pba.intervals import Interval
 from bisect import bisect_left
+import numpy as np
+import matplotlib.pyplot as plt
 
 
 class Joint(ABC):
@@ -119,6 +121,13 @@ def area_metric_sample(a: ArrayLike, b: ArrayLike):
     return sps.wasserstein_distance(a, b)
 
 
+#! not in use.
+def area_metric_np_numbers(a, b):
+    """when a and b are both numpy arrays of scalar numbers, compute the area metric accordingly"""
+    assert np.isscalar(a) and np.isscalar(b), "Both a and b must be scalar numbers."
+    return abs(a - b)
+
+
 def area_metric_number(a: Pbox | Number, b: Pbox | Number) -> float:
     """if any of a or b is a number, compute area metric accordingly"""
     from pyuncertainnumber import pba
@@ -226,10 +235,6 @@ def am_diff_register(A, B, debug=False):
     return result
 
 
-import numpy as np
-import matplotlib.pyplot as plt
-
-
 def intervals_from_res(res):
     """
     From a structured array `res` with fields: 'distance', 'x1', 'x2',
@@ -246,10 +251,6 @@ def intervals_from_res(res):
     return mask, intervals
 
 
-import numpy as np
-import matplotlib.pyplot as plt
-
-
 def plot_intervals_from_res(
     res, p, *, show_points=False, title=None, include_ties=False, ax=None
 ):
@@ -317,9 +318,6 @@ def plot_intervals_from_res(
     return ax
 
 
-import numpy as np
-
-
 def integrate_distance(res, p):
     """
     Integrate distance vs p with the trapezoidal rule,
@@ -531,3 +529,28 @@ def slide_pbox_towards_scalar(a, b):
             d2 = proposal_dd
 
     return Staircase(a.left - d1, a.right + d2)
+
+
+def double_metric(p: Number | Pbox | Interval, o: Number | Pbox | Interval):
+    """Double metric for two uncertain numbers.
+
+    args:
+        p: a prediction uncertain number (Pbox)
+        o: an observation uncertain number (Pbox or scalar)
+
+    note:
+        Typical case is for validation between prediction and observation where both are uncertain numbers.
+
+    """
+    if isinstance(p, Number):
+        p = Interval(p)
+
+    if isinstance(o, Number):
+        o = Interval(o)
+
+    return area_metric(p.left, o.left), area_metric(p.right, o.right)
+
+
+def conformal_double_metric(p: Number | Pbox | Interval, o: Number | Pbox | Interval):
+    """Propsed conformal version of the double metric, which takes the maximum of the two area metrics."""
+    return max(double_metric(p, o))
diff --git a/src/pyuncertainnumber/pba/intervals/intervalOperators.py b/src/pyuncertainnumber/pba/intervals/intervalOperators.py
index 8e053c32..185aa6dc 100644
--- a/src/pyuncertainnumber/pba/intervals/intervalOperators.py
+++ b/src/pyuncertainnumber/pba/intervals/intervalOperators.py
@@ -166,7 +166,8 @@ def make_vec_interval(vec):
     """
     from ...characterisation.uncertainNumber import UncertainNumber as UN
 
-    assert len(vec) > 1, "Interval must have more than one element"
+    # TODO: this is bugg. Interval vector has len(x)=1
+    # assert len(vec) > 1, "Interval must have more than one element"
 
     if isinstance(vec, Interval):
         return vec
diff --git a/src/pyuncertainnumber/pba/intervals/number.py b/src/pyuncertainnumber/pba/intervals/number.py
index 11526159..2de32197 100644
--- a/src/pyuncertainnumber/pba/intervals/number.py
+++ b/src/pyuncertainnumber/pba/intervals/number.py
@@ -146,7 +146,7 @@ def __getitem__(self, i: Union[int, slice]):  # make class indexable
     def to_numpy(self) -> np.ndarray:
         """transform interval objects to numpy arrays"""
         if self.scalar:
-            return np.array([self.lo.item(), self.hi.item()])
+            return np.array([self.lo, self.hi])
         else:
             return np.asarray((self.lo, self.hi)).T
 
@@ -222,14 +222,20 @@ def ravel(self):
 
     @property
     def lo(self) -> Union[ndarray, float]:
-        return self._lo
+        if self.scalar:
+            return self._lo.item()
+        else:
+            return self._lo
 
     # if len(self.shape)==0: return self._lo
     # return self._lo # return transpose(transpose(self.__val)[0]) # from shape (3,7,2) to (2,7,3) to (3,7)
 
     @property
     def hi(self) -> Union[ndarray, float]:
-        return self._hi
+        if self.scalar:
+            return self._hi.item()
+        else:
+            return self._hi
 
     @property
     def left(self):
@@ -397,9 +403,7 @@ def __rtruediv__(self, left):
         leftType = left.__class__.__name__
         # lo,hi = numpy.empty(self._lo.shape),numpy.empty(self._hi.shape)
         self_lo, self_hi = self.lo, self.hi
-        self_straddle_zero = numpy.any(
-            (self_lo.flatten() <= 0) & (self_hi.flatten() >= 0)
-        )
+        self_straddle_zero = numpy.any((self_lo <= 0) & (self_hi >= 0))
         if self_straddle_zero:
             raise ZeroDivisionError
         if (leftType == "ndarray") | (leftType in NUMERIC_TYPES):
@@ -467,6 +471,33 @@ def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
         if "out" in kwargs and kwargs["out"] is not None:
             return NotImplemented
 
+        binary = {
+            np.add: ("__add__", "__radd__"),
+            np.subtract: ("__sub__", "__rsub__"),
+            np.multiply: ("__mul__", "__rmul__"),
+            np.true_divide: ("__truediv__", "__rtruediv__"),
+            np.floor_divide: ("__floordiv__", "__rfloordiv__"),
+            np.power: ("__pow__", "__rpow__"),
+            np.maximum: ("__max__", "__rmax__"),
+            np.minimum: ("__min__", "__rmin__"),
+        }
+
+        if ufunc in binary and len(inputs) == 2:
+            left, right = inputs
+            l_name, r_name = binary[ufunc]
+
+            if left is self:
+                # self (op) right
+                return getattr(self, l_name)(right)
+            elif right is self:
+                # left (op) self
+                return getattr(self, r_name)(left)
+            else:
+                return NotImplemented
+
+        if len(inputs) != 1 or inputs[0] is not self:
+            return NotImplemented
+
         if ufunc is np.sin:
             return self.sin()
         if ufunc is np.cos:
diff --git a/src/pyuncertainnumber/pba/pbox_abc.py b/src/pyuncertainnumber/pba/pbox_abc.py
index 18224e96..539bfdb9 100644
--- a/src/pyuncertainnumber/pba/pbox_abc.py
+++ b/src/pyuncertainnumber/pba/pbox_abc.py
@@ -863,7 +863,9 @@ def display(self, *args, **kwargs):
         self.plot(*args, **kwargs)
         plt.show()
 
-    def plot_probability_bound(self, x: float, ax=None, **kwargs):
+    def plot_probability_bound(
+        self, x: float, ax=None, linecolor="r", markercolor="r", **kwargs
+    ):
         """plot the probability bound at a certain quantile x
 
         note:
@@ -880,12 +882,12 @@ def plot_probability_bound(self, x: float, ax=None, **kwargs):
         ax.plot(
             [x, x],
             [p_lo, p_hi],
-            c="r",
+            c=linecolor,
             label="probability bound",
             zorder=50,
         )
-        ax.scatter(x, p_lo, c="r", marker="^", zorder=50)
-        ax.scatter(x, p_hi, c="r", marker="v", zorder=50)
+        ax.scatter(x, p_lo, c=markercolor, marker="^", zorder=50)
+        ax.scatter(x, p_hi, c=markercolor, marker="v", zorder=50)
         return ax
 
     def plot_quantile_bound(self, p: float, ax=None, **kwargs):
@@ -982,9 +984,34 @@ def __rpow__(self, other: Number):
     def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
         if method != "__call__":
             return NotImplemented
-        if len(inputs) != 1 or inputs[0] is not self:
+        if kwargs.get("out", None) is not None:
             return NotImplemented
-        if "out" in kwargs and kwargs["out"] is not None:
+
+        binary = {
+            np.add: ("__add__", "__radd__"),
+            np.subtract: ("__sub__", "__rsub__"),
+            np.multiply: ("__mul__", "__rmul__"),
+            np.true_divide: ("__truediv__", "__rtruediv__"),
+            np.floor_divide: ("__floordiv__", "__rfloordiv__"),
+            np.power: ("__pow__", "__rpow__"),
+            np.maximum: ("__max__", "__rmax__"),
+            np.minimum: ("__min__", "__rmin__"),
+        }
+
+        if ufunc in binary and len(inputs) == 2:
+            left, right = inputs
+            l_name, r_name = binary[ufunc]
+
+            if left is self:
+                # self (op) right
+                return getattr(self, l_name)(right)
+            elif right is self:
+                # left (op) self
+                return getattr(self, r_name)(left)
+            else:
+                return NotImplemented
+
+        if len(inputs) != 1 or inputs[0] is not self:
             return NotImplemented
 
         if ufunc is np.sin:
@@ -1035,6 +1062,26 @@ def sample(self, n_sam):
         alpha = np.squeeze(qmc.LatinHypercube(d=1).random(n=n_sam))
         return self.alpha_cut(alpha)
 
+    def precise_sample(
+        self,
+        n_a: int,
+        theta: float = None,
+        n_e: int = None,
+    ):
+        """Generate precise samples from a p-box"""
+        if (theta is None) and (n_e is None):
+            raise ValueError("Either theta or n_e must be provided.")
+        if theta is not None and n_e is None:
+            assert 0 <= theta <= 1, "Theta must be in the range [0, 1]."
+            focal_elements = self.sample(n_a)
+            return (focal_elements.hi - focal_elements.lo) * theta + focal_elements.lo
+        if n_e is not None and theta is None:
+            theta = np.random.uniform(0, 1, size=n_e)
+            focal_elements = self.sample(n_a)
+            return (focal_elements.hi - focal_elements.lo)[None, :] * theta[
+                :, None
+            ] + focal_elements.lo[None, :]
+
     def discretise(self, n=None) -> Interval:
         """alpha-cut discretisation of the p-box without outward rounding
 
diff --git a/src/pyuncertainnumber/propagation/taylor_expansion.py b/src/pyuncertainnumber/propagation/taylor_expansion.py
index 5aceb1bb..6f420e0a 100644
--- a/src/pyuncertainnumber/propagation/taylor_expansion.py
+++ b/src/pyuncertainnumber/propagation/taylor_expansion.py
@@ -1,10 +1,5 @@
 import os
 
-os.environ["JAX_PLATFORMS"] = "cpu"
-import jax
-import jax.numpy as jnp
-
-
 """ Taylor expansions for the moments of functions of random variables """
 
 
@@ -36,6 +31,8 @@ def taylor_expansion_method(func, mean, *, var=None, cov=None) -> tuple:
         >>> mu_, var_ = taylor_expansion_method(func=bar, mean=MEAN, cov=COV)
 
     """
+    os.environ["JAX_PLATFORMS"] = "cpu"
+
     if mean.ndim == 1:  # random vector
         return taylor_expansion_method_vector(func, mean, cov)
     elif mean.ndim == 0:  # scalar random variable
@@ -45,6 +42,9 @@ def taylor_expansion_method(func, mean, *, var=None, cov=None) -> tuple:
 def taylor_expansion_method_scalar(func, mean, var) -> tuple:
     """For scalar random variable only"""
 
+    import jax
+    import jax.numpy as jnp
+
     # gradient
     d1f = jax.grad(func)(mean)
 
@@ -61,6 +61,9 @@ def taylor_expansion_method_scalar(func, mean, var) -> tuple:
 def taylor_expansion_method_vector(func, mean, cov) -> tuple:
     """For random vector only"""
 
+    import jax
+    import jax.numpy as jnp
+
     # gradient
     d1f = jax.grad(func)(mean)
 
diff --git a/tests/test_double_metric.py b/tests/test_double_metric.py
new file mode 100644
index 00000000..ae4db2e2
--- /dev/null
+++ b/tests/test_double_metric.py
@@ -0,0 +1,38 @@
+from pyuncertainnumber import pba
+from pyuncertainnumber.pba.core import double_metric, conformal_double_metric
+
+
+def test_double_metric():
+    x_left = 7
+    y_left = pba.I(14, 19)
+
+    x_middle = pba.I(4, 9)
+    y_middle = pba.I(13, 18)
+
+    x_right = pba.I(3, 11)
+    y_right = pba.I(8, 17)
+
+    for a, b in [(x_left, y_left), (x_middle, y_middle), (x_right, y_right)]:
+        print(double_metric(a, b))
+
+    assert double_metric(x_left, y_left) == (7, 12)
+    assert double_metric(x_middle, y_middle) == (9, 9)
+    assert double_metric(x_right, y_right) == (5, 6)
+
+
+def test_conformal_double_metric():
+    x_left = 7
+    y_left = pba.I(14, 19)
+
+    x_middle = pba.I(4, 9)
+    y_middle = pba.I(13, 18)
+
+    x_right = pba.I(3, 11)
+    y_right = pba.I(8, 17)
+
+    for a, b in [(x_left, y_left), (x_middle, y_middle), (x_right, y_right)]:
+        print("conformal version", conformal_double_metric(a, b))
+
+    assert conformal_double_metric(x_left, y_left) == 12
+    assert conformal_double_metric(x_middle, y_middle) == 9
+    assert conformal_double_metric(x_right, y_right) == 6
diff --git a/tests/test_numpy_binary_compatability.py b/tests/test_numpy_binary_compatability.py
new file mode 100644
index 00000000..e69de29b