unumpy.average: init

Ref: lmfit#38

CHANGES.rst

@@ -7,6 +7,7 @@ Unreleased
 Fixes:
 
 - fix `readthedocs` configuration so that the build passes (#254)
+- Add ``unumpy.average`` to calculate uncertainties aware average (#38)
 
 3.2.2   2024-July-08
 -----------------------

doc/numpy_guide.rst

@@ -144,6 +144,22 @@ functions is available in the documentation for  :mod:`uncertainties.umath`.
 .. index::
    pair: testing and operations (in arrays); NaN
 
+Uncertainties aware average
+---------------------------
+
+If you have measured a certain value multiple times, with a different uncertainty every measurement. Averaging over the results in a manner aware of the different uncertainties, is not trivial. The function ``unumpy.average()`` does that:
+
+>>> measurements = numpy.array([2.1, 2.0, 2.05, 2.08, 2.02])
+>>> stds = numpy.array([0.05, 0.03, 0.04, 0.06, 0.05])
+>>> arr = unumpy.uarray(measurements, stds)
+>>> unumpy.average(arr)
+2.03606+/-0.00018
+
+Note how that function gives a value different from numpy's ``mean`` function:
+
+>>> numpy.mean(arr)
+2.050+/-0.021
+
 NaN testing and NaN-aware operations
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 

tests/test_unumpy.py

@@ -1,3 +1,5 @@
+import pytest
+
 try:
     import numpy
 except ImportError:
@@ -300,3 +302,9 @@ def test_array_comparisons():
     # For matrices, 1D arrays are converted to 2D arrays:
     mat = unumpy.umatrix([1, 2], [1, 4])
     assert numpy.all((mat == [mat[0, 0], 4]) == [True, False])
+
+
+def test_average_type_check():
+    arr = numpy.array(["bla"])
+    with pytest.raises(ValueError):
+        unumpy.average(arr)

uncertainties/unumpy/core.py

@@ -30,10 +30,30 @@
     # Utilities:
     "nominal_values",
     "std_devs",
+    "average",
     # Classes:
     "matrix",
 ]
 
+
+def average(arr, axis=None):
+    """
+    Return a weighted averaged along with a weighted mean over a certain axis.
+
+    By default, operates on all axes of the given array.
+    """
+    if arr.dtype != numpy.dtype("object") or (
+        not isinstance(arr.flat[0], uncert_core.Variable)
+    ):
+        raise ValueError(
+            "unumpy.average is meant to operate only upon numpy arrays of ufloats."
+        )
+    weights = std_devs(arr) ** -2
+    weights_sum = weights.sum(axis=axis)
+    weighted_mean = (nominal_values(arr) * weights).sum(axis=axis) / weights_sum
+    return uarray(weighted_mean, weights_sum**-1 / 2)
+
+
 ###############################################################################
 # Utilities: