documentation

Author: renecotyfanboy

docs/examples/plot_5_emcee_arviz_numpyro.py

@@ -1,5 +1,5 @@
 """
-# Using external samples easily
+# Using external samples
 
 `emcee`, `arviz`, and `numpyro` are all popular MCMC packages. ChainConsumer
 provides class methods to turn results from these packages into chains efficiently.

docs/examples/plot_7_multimodal_chains.py

@@ -0,0 +1,84 @@
+"""
+# Multimodal distributions
+
+`emcee`, `arviz`, and `numpyro` are all popular MCMC packages. ChainConsumer
+provides class methods to turn results from these packages into chains efficiently.
+
+If you want to request support for another type of chain, please open a
+[discussion](https://github.com/Samreay/ChainConsumer/discussions) with a code
+example, and we can add it in. The brave may even provide a PR!
+"""
+
+import numpy as np
+import pandas as pd
+
+from chainconsumer import Chain, ChainConsumer
+from chainconsumer.statistics import SummaryStatistic
+
+# %%
+# First, let's build some dummy data
+
+rng = np.random.default_rng(42)
+size = 60_000
+
+eta = rng.normal(loc=0.0, scale=0.8, size=size)
+
+phi = np.asarray(
+    [rng.gamma(shape=2.5, scale=0.4, size=size // 2) - 3.0, 3.0 - rng.gamma(shape=5.0, scale=0.35, size=(size // 2))]
+).flatten()
+
+rng.shuffle(phi)
+
+df = pd.DataFrame({"eta": eta, "phi": phi})
+
+# %%
+# To build a multimodal chain, you simply have to pass `multimodal=True` when building the chain. To work, it requires
+# you to specify `SummaryStatistic.HDI` as the summary statistic.
+
+chain_multimodal = Chain(
+    samples=df.copy(),
+    name="posterior-multimodal",
+    statistics=SummaryStatistic.HDI,
+    multimodal=True,  # <- Here
+)
+
+# %%
+# Now, if you add this `Chain` to a plotter, it will try to look for sub-intervals and display them.
+
+cc = ChainConsumer()
+cc.add_chain(chain_multimodal)
+fig = cc.plotter.plot()
+
+# %%
+# Let's compare with what would happen if you don't use a multimodal chain. We use the same data as before but don't
+# warn `ChainConsumer` that we expect the chains to be multimodal.
+
+chain_unimodal = Chain(samples=df.copy(), name="posterior-unimodal", statistics=SummaryStatistic.HDI, multimodal=False)
+
+cc.add_chain(chain_unimodal)
+fig = cc.plotter.plot()
+
+# %%
+# Let's try with even more modes.
+
+eta = np.asarray(
+    [
+        rng.normal(loc=-3, scale=0.8, size=size // 3),
+        rng.normal(loc=0.0, scale=0.8, size=size // 3),
+        rng.normal(loc=+3, scale=0.8, size=size // 3),
+    ]
+).flatten()
+
+
+rng.shuffle(eta)
+
+df = pd.DataFrame({"eta": eta, "phi": phi})
+
+chain_multimodal = Chain(
+    samples=df.copy(), name="posterior-multimodal", statistics=SummaryStatistic.HDI, multimodal=True
+)
+
+cc = ChainConsumer()
+cc.add_chain(chain_multimodal)
+fig = cc.plotter.plot()
+fig.tight_layout()

src/chainconsumer/analysis.py

@@ -41,7 +41,7 @@ def _mask_to_intervals(
     if mask[-1]:
         ends = np.concatenate((ends, len(mask) - 1))
 
-    intervals = [(float(x[s]), float(x[e])) for s, e in zip(starts, ends, strict=False) if x[e] > x[s]]
+    intervals = [(float(x[s]), float(x[e])) for s, e in zip(starts, ends, strict=True) if x[e] > x[s]]
 
     return intervals
 
@@ -202,8 +202,8 @@ def get_summary(
 
         Returns:
             Mapping from chain name to marginal summaries. Each entry is either a single
-            :class:`~chainconsumer.analysis.Bound` for unimodal chains or a list of bounds
-            describing the disjoint HDI bands when ``chain.multimodal`` is set.
+            [`Bound`][chainconsumer.analysis.Bound] for unimodal chains or a list of bounds
+            describing the disjoint HDI bands when `chain.multimodal` is set.
         """
         results = {}
         if chains is None:

src/chainconsumer/plotter.py

@@ -961,59 +961,60 @@ def _plot_bars(
         if chain.bar_shade:
             base_bound = self.parent.analysis.get_parameter_summary(chain, column)
 
-            if chain.multimodal:
-                intervals = self.parent.analysis.get_parameter_hdi_intervals(chain, column)
-                display_bounds = self.parent.analysis.get_parameter_multimodal_bounds(
-                    chain,
-                    column,
-                    intervals,
-                )
-
-                # If we get a single interval, fallback to unimodal HDI
-                if len(intervals) < 2:
-                    intervals = [(base_bound.lower, base_bound.upper)]
-
-            else:
-                display_bounds = base_bound
-                intervals = [(display_bounds.lower, display_bounds.upper)]
-                intervals = np.clip(intervals, a_min=xs.min(), a_max=xs.max())
-
-            for lower, upper_ in intervals:
-                x = np.linspace(lower, upper_, 1000)
-
-                if flip:
-                    ax.fill_betweenx(
-                        x,
-                        np.zeros_like(x),
-                        interpolator(x),
-                        color=chain.color,
-                        alpha=0.2,
-                        zorder=chain.zorder,
+            if base_bound is not None and base_bound.lower is not None and base_bound.upper is not None:
+                if chain.multimodal:
+                    intervals = self.parent.analysis.get_parameter_hdi_intervals(chain, column)
+                    display_bounds = self.parent.analysis.get_parameter_multimodal_bounds(
+                        chain,
+                        column,
+                        intervals,
                     )
 
+                    # If we get a single interval, fallback to unimodal HDI
+                    if len(intervals) < 2:
+                        intervals = [(base_bound.lower, base_bound.upper)]
+
                 else:
-                    ax.fill_between(
-                        x,
-                        np.zeros_like(x),
-                        interpolator(x),
-                        color=chain.color,
-                        alpha=0.2,
-                        zorder=chain.zorder,
-                    )
+                    display_bounds = base_bound
+                    intervals = [(display_bounds.lower, display_bounds.upper)]
+                    intervals = np.clip(intervals, a_min=xs.min(), a_max=xs.max())
+
+                for lower, upper_ in intervals:
+                    x = np.linspace(lower, upper_, 1000)
+
+                    if flip:
+                        ax.fill_betweenx(
+                            x,
+                            np.zeros_like(x),
+                            interpolator(x),
+                            color=chain.color,
+                            alpha=0.2,
+                            zorder=chain.zorder,
+                        )
 
-            if summary:
-                label = self.config.get_label(column)
-                label_core = label.strip("$") if isinstance(column, str) else None
-                label_text = label_core or None
+                    else:
+                        ax.fill_between(
+                            x,
+                            np.zeros_like(x),
+                            interpolator(x),
+                            color=chain.color,
+                            alpha=0.2,
+                            zorder=chain.zorder,
+                        )
 
-                title = self.parent.analysis.get_parameter_text(
-                    display_bounds,
-                    wrap=True,
-                    label=label_text,
-                )
+                if summary:
+                    label = self.config.get_label(column)
+                    label_core = label.strip("$") if isinstance(column, str) else None
+                    label_text = label_core or None
+
+                    title = self.parent.analysis.get_parameter_text(
+                        display_bounds,
+                        wrap=True,
+                        label=label_text,
+                    )
 
-                if title:
-                    ax.set_title(title, fontsize=self.config.summary_font_size)
+                    if title:
+                        ax.set_title(title, fontsize=self.config.summary_font_size)
 
         return float(ys.max())