ENH: add Moran.plot_simulation and Moran_BV.plot_simulation

ci/310-numba-oldest.yaml

@@ -20,6 +20,7 @@ dependencies:
   - folium
   - mapclassify
   - matplotlib
+  - seaborn
   - pytest
   - pytest-cov
   - pytest-xdist

ci/310-oldest.yaml

@@ -19,6 +19,7 @@ dependencies:
   - folium
   - mapclassify
   - matplotlib
+  - seaborn
   - pytest
   - pytest-cov
   - pytest-xdist

ci/311-latest.yaml

@@ -15,6 +15,7 @@ dependencies:
   - folium
   - mapclassify
   - matplotlib
+  - seaborn
   - pytest
   - pytest-cov
   - pytest-xdist

ci/311-numba-latest.yaml

@@ -16,6 +16,7 @@ dependencies:
   - folium
   - mapclassify
   - matplotlib
+  - seaborn
   - pytest
   - pytest-cov
   - pytest-xdist

ci/312-dev.yaml

@@ -12,6 +12,7 @@ dependencies:
   - folium
   - mapclassify
   - matplotlib
+  - seaborn
   - pytest
   - pytest-cov
   - pytest-xdist

ci/312-latest.yaml

@@ -16,6 +16,7 @@ dependencies:
   - folium
   - mapclassify
   - matplotlib
+  - seaborn
   - pytest
   - pytest-cov
   - pytest-xdist

ci/312-numba-dev.yaml

@@ -13,6 +13,7 @@ dependencies:
   - folium
   - mapclassify
   - matplotlib
+  - seaborn
   - pytest
   - pytest-cov
   - pytest-xdist

ci/312-numba-latest.yaml

@@ -16,6 +16,7 @@ dependencies:
   - folium
   - mapclassify
   - matplotlib
+  - seaborn
   - pytest
   - pytest-cov
   - pytest-xdist

esda/moran.py

@@ -352,6 +352,58 @@ def plot_scatter(
             fitline_kwds=fitline_kwds,
         )
 
+    def plot_simulation(self, ax=None, legend=False, fitline_kwds=None, **kwargs):
+        """
+        Global Moran's I simulated reference distribution.
+
+        Parameters
+        ----------
+        ax : matplotlib.axes.Axes, optional
+            Pre-existing axes for the plot, by default None.
+        legend : bool, optional
+            Plot a legend, by default False
+        fitline_kwds : dict, optional
+            Additional keyword arguments for vertical Moran fit line, by default None.
+        **kwargs : keyword arguments, optional
+            Additional keyword arguments for KDE plot passed to ``seaborn.kdeplot``,
+            by default None.
+
+        Returns
+        -------
+        matplotlib.axes.Axes
+            Axes object with the Moran scatterplot.
+
+        Notes
+        -----
+        This requires optional dependencies ``matplotlib`` and ``seaborn``.
+
+        Examples
+        --------
+        >>> import libpysal
+        >>> w = libpysal.io.open(libpysal.examples.get_path("stl.gal")).read()
+        >>> f = libpysal.io.open(libpysal.examples.get_path("stl_hom.txt"))
+        >>> y = np.array(f.by_col['HR8893'])
+        >>> from esda.moran import Moran
+        >>> mi = Moran(y,  w)
+
+        Default plot:
+
+        >>> mi.plot_simulation()
+
+        Customized styling that turns the distribution into a pink line and line
+        indicating I to a black line:
+
+        >>> mi.plot_simulation(fitline_kwds={"color": "k"}, color="pink", shade=False)
+        """
+        return _simulation_plot(
+            self,
+            ax=ax,
+            legend=legend,
+            bivariate=False,
+            fitline_kwds=fitline_kwds,
+            **kwargs,
+        )
+
 
 class Moran_BV:  # noqa: N801
     """
@@ -565,6 +617,40 @@ def by_col(
             **stat_kws,
         )
 
+    def plot_simulation(self, ax=None, legend=False, fitline_kwds=None, **kwargs):
+        """
+        Global Moran's I simulated reference distribution.
+
+        Parameters
+        ----------
+        ax : matplotlib.axes.Axes, optional
+            Pre-existing axes for the plot, by default None.
+        legend : bool, optional
+            Plot a legend, by default False
+        fitline_kwds : dict, optional
+            Additional keyword arguments for vertical Moran fit line, by default None.
+        **kwargs : keyword arguments, optional
+            Additional keyword arguments for KDE plot passed to ``seaborn.kdeplot``,
+            by default None.
+
+        Returns
+        -------
+        matplotlib.axes.Axes
+            Axes object with the Moran scatterplot.
+
+        Notes
+        -----
+        This requires optional dependencies ``matplotlib`` and ``seaborn``.
+        """
+        return _simulation_plot(
+            self,
+            ax=ax,
+            legend=legend,
+            bivariate=True,
+            fitline_kwds=fitline_kwds,
+            **kwargs,
+        )
+
 
 def Moran_BV_matrix(variables, w, permutations=0, varnames=None):  # noqa: N802
     """
@@ -2024,6 +2110,49 @@ def _scatterplot(
     return ax
 
 
+def _simulation_plot(
+    moran, ax=None, legend=False, bivariate=False, fitline_kwds=None, **kwargs
+):
+    try:
+        import seaborn as sns
+        from matplotlib import pyplot as plt
+    except ImportError as err:
+        raise ImportError(
+            "matplotlib and seaborn must be installed to plot the simulation."
+        ) from err
+    # to set default as an empty dictionary that is later filled with defaults
+    if fitline_kwds is None:
+        fitline_kwds = dict()
+
+    if ax is None:
+        _, ax = plt.subplots()
+
+    # plot distribution
+    shade = kwargs.pop("shade", True)
+    color = kwargs.pop("color", "#bababa")
+    sns.kdeplot(
+        moran.sim,
+        fill=shade,
+        color=color,
+        ax=ax,
+        label="Distribution of simulated Is",
+        **kwargs,
+    )
+
+    exp = moran.EI_sim if bivariate else moran.EI
+
+    # customize plot
+    fitline_kwds.setdefault("color", "#d6604d")
+    ax.vlines(moran.I, 0, 1, **fitline_kwds, label="Moran's I")
+    ax.vlines(exp, 0, 1, label="Expected I")
+    ax.set_title("Reference Distribution")
+    ax.set_xlabel(f"Moran's I: {moran.I:.2f}")
+
+    if legend:
+        ax.legend()
+    return ax
+
+
 # --------------------------------------------------------------
 # Conditional Randomization Moment Estimators
 # --------------------------------------------------------------

esda/tests/test_moran.py

@@ -115,6 +115,95 @@ def test_by_col(self):
         np.testing.assert_allclose(sidr, 0.24772519320480135, atol=ATOL, rtol=RTOL)
         np.testing.assert_allclose(pval, 0.001)
 
+    @parametrize_sac
+    def test_plot_simulation(self, w):
+        pytest.importorskip("seaborn")
+
+        m = moran.Moran(sac1.WHITE, w=w)
+        ax = m.plot_simulation()
+
+        assert len(ax.collections) == 3
+
+        kde = ax.collections[0]
+        np.testing.assert_array_almost_equal(
+            kde.get_facecolor(),
+            [[0.7294117647058823, 0.7294117647058823, 0.7294117647058823, 0.25]],
+        )
+        assert kde.get_fill()
+        assert len(kde.get_paths()[0]) == 403
+
+        i_vline = ax.collections[1]
+        np.testing.assert_array_almost_equal(
+            i_vline.get_color(),
+            [[0.8392156862745098, 0.3764705882352941, 0.30196078431372547, 1.0]],
+        )
+        assert i_vline.get_label() == "Moran's I"
+        np.testing.assert_array_almost_equal(
+            i_vline.get_paths()[0].vertices,
+            np.array([[m.I, 0.0], [m.I, 1.0]]),
+        )
+
+        ei_vline = ax.collections[2]
+        np.testing.assert_array_almost_equal(
+            ei_vline.get_color(),
+            [[0.12156863, 0.46666667, 0.70588235, 1.0]],
+        )
+        assert ei_vline.get_label() == "Expected I"
+        np.testing.assert_array_almost_equal(
+            ei_vline.get_paths()[0].vertices,
+            np.array([[m.EI, 0.0], [m.EI, 1.0]]),
+        )
+
+    @parametrize_sac
+    def test_plot_simulation_custom(self, w):
+        pytest.importorskip("seaborn")
+        plt = pytest.importorskip("matplotlib.pyplot")
+
+        m = moran.Moran(sac1.WHITE, w=w)
+
+        _, ax = plt.subplots(figsize=(12, 12))
+        ax = m.plot_simulation(
+            ax=ax, fitline_kwds={"color": "red"}, color="pink", shade=False, legend=True
+        )
+
+        assert len(ax.collections) == 2
+        assert len(ax.lines) == 1
+
+        kde = ax.lines[0]
+        np.testing.assert_array_almost_equal(
+            kde.get_color(),
+            [1.0, 0.75294118, 0.79607843, 1],
+        )
+        assert len(kde.get_path()) == 200
+
+        i_vline = ax.collections[0]
+        np.testing.assert_array_almost_equal(
+            i_vline.get_color(),
+            [[1.0, 0.0, 0.0, 1.0]],
+        )
+        assert i_vline.get_label() == "Moran's I"
+        np.testing.assert_array_almost_equal(
+            i_vline.get_paths()[0].vertices,
+            np.array([[m.I, 0.0], [m.I, 1.0]]),
+        )
+
+        ei_vline = ax.collections[1]
+        np.testing.assert_array_almost_equal(
+            ei_vline.get_color(),
+            [[0.12156863, 0.46666667, 0.70588235, 1.0]],
+        )
+        assert ei_vline.get_label() == "Expected I"
+        np.testing.assert_array_almost_equal(
+            ei_vline.get_paths()[0].vertices,
+            np.array([[m.EI, 0.0], [m.EI, 1.0]]),
+        )
+
+        assert ax.get_legend_handles_labels()[1] == [
+            "Distribution of simulated Is",
+            "Moran's I",
+            "Expected I",
+        ]
+
     @parametrize_sac
     def test_Moran_plot_scatter(self, w):
         import matplotlib

pyproject.toml

@@ -49,6 +49,7 @@ plus = [
     "matplotlib",
     "numba>=0.57",
     "rtree>=1.0",
+    "seaborn",
 ]
 tests = [
     "codecov",
@@ -72,7 +73,6 @@ docs = [
 ]
 notebooks = [
     "matplotlib-scalebar",
-    "seaborn",
     "watermark",
 ]