Add from_stack

docs/create.rst

@@ -78,6 +78,13 @@ From publication
 TODO: show example how to take covar (or par errors) from a
 publication or blog post, i.e. as inputs.
 
+.. _create_from_stack:
+
+From stack
+----------
+
+TODO: document :meth:`MultiNorm.from_stack`
+
 .. _create_from_product:
 
 From product

docs/theory.rst

@@ -166,6 +166,13 @@ Add note that for MVN the covar matrix for conditional doesn't depend on paramet
 
 TODO: document and make example in the analyse section using iminuit.
 
+.. _theory_stack:
+
+Stacked distribution
+--------------------
+
+TODO: document :meth:`MultiNorm.from_stack`
+
 .. _theory_product:
 
 Product distribution

multinorm.py

@@ -12,7 +12,7 @@
 from pkg_resources import get_distribution, DistributionNotFound
 import numpy as np
 import pandas as pd
-from scipy.linalg import eigh
+from scipy.linalg import eigh, block_diag
 from scipy.stats import multivariate_normal
 
 __all__ = ["MultiNorm"]
@@ -142,14 +142,43 @@ def from_samples(cls, points, names=None):
         cov = np.cov(points, rowvar=False)
         return cls(mean, cov, names)
 
+    @classmethod
+    def from_stack(cls, distributions):
+        """Create `MultiNorm` as stacked distributions.
+
+        Stacking means the ``names`` and ``mean`` vectors
+        are concatenated, and the ``cov`` matrices are
+        combined into a block diagonal matrix, with zeros
+        for the off-diagonal parts.
+
+        This represents the combined measurement, assuming
+        the individual distributions are for different parameters.
+
+        See :ref:`create_from_stack` and :ref:`theory_stack`.
+
+        Parameters
+        ----------
+        distributions : list
+            Python list of `MultiNorm` distributions.
+
+        Returns
+        -------
+        MultiNorm
+            Stacked distribution
+        """
+        names = np.concatenate([_.names for _ in distributions])
+        cov = block_diag(*[_.cov for _ in distributions])
+        mean = np.concatenate([_.mean for _ in distributions])
+        return cls(mean, cov, names)
+
     @classmethod
     def from_product(cls, distributions):
         """Create `MultiNorm` as product distribution.
 
         This represents the joint likelihood distribution, assuming
         the individual distributions are from independent measurements.
 
-        See :ref:`theory_product` .
+        See :ref:`create_from_product` and :ref:`theory_product` .
 
         Parameters
         ----------
@@ -561,7 +590,8 @@ def __init__(self, names, n):
             if len(names) != n:
                 raise ValueError(f"len(names) = {len(names)}, expected n={n}")
 
-        self.names = names
+        # TODO: change to Numpy array!
+        self.names = list(names)
 
     def get_idx(self, pars):
         """Create parameter index array.

test_multinorm.py

@@ -129,6 +129,17 @@ def test_from_samples():
     assert_allclose(mn.mean, [11, 20, 30])
     assert_allclose(mn.cov, [[2, 0, 0], [0, 0, 0], [0, 0, 0]])
 
+def test_from_stack():
+    d1 = MultiNorm(mean=[1, 2], cov=np.full((2, 2), 2), names=["a", "b"])
+    d2 = MultiNorm(mean=[3, 4, 5], cov=np.full((3, 3), 3), names=["c", "d", "e"])
+
+    mn = MultiNorm.from_stack([d1, d2])
+
+    assert mn.names == ["a", "b", "c", "d", "e"]
+    assert_allclose(mn.mean, [1, 2, 3, 4, 5])
+    assert_allclose(mn.cov.values[0], [2, 2, 0, 0, 0])
+    assert_allclose(mn.cov.values[4], [0, 0, 3, 3, 3])
+
 
 def test_from_product():
     d1 = MultiNorm(mean=[0, 0], names=["a", "b"])