Repository: arviz-devs/arviz

Issue:
Add plot_ecdf
## Tell us about it

It would be nice to have an empirical CDF (ECDF) plot function that isn't tied to `plot_ppc`. As far as information a `plot_ecdf` might display, I'm thinking it would behave pretty similarly to `plot_density`, but with an ECDF. One extra thing you might want to do is specify some quantiles that would be drawn on the plot. The most common style for ECDF is the line visualization.

![ecdf](https://user-images.githubusercontent.com/8673634/83200946-4ac5cf80-a0f9-11ea-8f96-dbd5e5c2a88b.png)

Another style that we might want to support is one that shows the intervals (e.g. see below from https://statisticsglobe.com/empirical-cumulative-distribution-function-ecdf-r)
![R-Basics-the-dim-Function-in-R](https://user-images.githubusercontent.com/8673634/83200030-41d3fe80-a0f7-11ea-9882-48f6877dbc79.png)
However, in my experience, this is only useful for small sample sizes, and in general for posterior samples, the line is cleaner.

## Thoughts on implementation

It's already possible to compute the ECDF (we already have `arviz.plots.backends.matplotlib.ppcplot._empirical_cdf` and
`arviz.plots.backends.bokeh.ppcplot._empirical_cdf`), so this should be pretty straightforward to implement.

Discussion:
Rishabh261998: Thanks, @OriolAbril  for the suggestions. First I am thinking to add all the computation code in the general function and writing the broken variation. After this, we can look into the performance of the compute_ecdf and also how numba can be used to speed up if time permits. Also, there is a more optimized code for computing gamma as presented in that paper, could also implement that.
Rishabh261998: @sethaxen  and @aloctavodia please go through it whenever you get time. Thanks
Rishabh261998: Thanks @sethaxen I have started making the suggested changes and will soon make a commit. 
codecov[bot]: # [Codecov](https://codecov.io/gh/arviz-devs/arviz/pull/1753?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=arviz-devs) Report
> Merging [#1753](https://codecov.io/gh/arviz-devs/arviz/pull/1753?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=arviz-devs) (de5fdd1) into [main](https://codecov.io/gh/arviz-devs/arviz/commit/928f8f4c34ae20f64bc787486f34dc0602bfae6f?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=arviz-devs) (928f8f4) will **decrease** coverage by `5.01%`.
> The diff coverage is `5.35%`.

[![Impacted file tree graph](https://codecov.io/gh/arviz-devs/arviz/pull/1753/graphs/tree.svg?width=650&height=150&src=pr&token=KZTdQCEf7t&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=arviz-devs)](https://codecov.io/gh/arviz-devs/arviz/pull/1753?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=arviz-devs)

```diff
@@            Coverage Diff             @@
##             main    #1753      +/-   ##
==========================================
- Coverage   91.18%   86.16%   -5.02%     
==========================================
  Files         115      118       +3     
  Lines       12248    12416     +168     
==========================================
- Hits        11168    10698     -470     
- Misses       1080     1718     +638     
```


| [Impacted Files](https://codecov.io/gh/arviz-devs/arviz/pull/1753?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=arviz-devs) | Coverage Δ | |
|---|---|---|
| [arviz/plots/backends/bokeh/ecdfplot.py](https://codecov.io/gh/arviz-devs/arviz/pull/1753/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=arviz-devs#diff-YXJ2aXovcGxvdHMvYmFja2VuZHMvYm9rZWgvZWNkZnBsb3QucHk=) | `0.00% <0.00%> (ø)` | |
| [arviz/plots/backends/matplotlib/ecdfplot.py](https://codecov.io/gh/arviz-devs/arviz/pull/1753/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=arviz-devs#diff-YXJ2aXovcGxvdHMvYmFja2VuZHMvbWF0cGxvdGxpYi9lY2RmcGxvdC5weQ==) | `0.00% <0.00%> (ø)` | |
| [arviz/plots/ecdfplot.py](https://codecov.io/gh/arviz-devs/arviz/pull/1753/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=arviz-devs#diff-YXJ2aXovcGxvdHMvZWNkZnBsb3QucHk=) | `7.92% <7.92%> (ø)` | |
| [arviz/plots/\_\_init\_\_.py](https://codecov.io/gh/arviz-devs/arviz/pull/1753/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=arviz-devs#diff-YXJ2aXovcGxvdHMvX19pbml0X18ucHk=) | `100.00% <100.00%> (ø)` | |
| [arviz/data/io\_pymc3\_3x.py](https://codecov.io/gh/arviz-devs/arviz/pull/1753/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=arviz-devs#diff-YXJ2aXovZGF0YS9pb19weW1jM18zeC5weQ==) | `18.18% <0.00%> (-72.06%)` | :arrow_down: |
| [arviz/data/io\_pystan.py](https://codecov.io/gh/arviz-devs/arviz/pull/1753/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=arviz-devs#diff-YXJ2aXovZGF0YS9pb19weXN0YW4ucHk=) | `62.12% <0.00%> (-33.96%)` | :arrow_down: |
| [arviz/data/io\_cmdstanpy.py](https://codecov.io/gh/arviz-devs/arviz/pull/1753/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=arviz-devs#diff-YXJ2aXovZGF0YS9pb19jbWRzdGFucHkucHk=) | `70.68% <0.00%> (-24.76%)` | :arrow_down: |
| [arviz/data/io\_emcee.py](https://codecov.io/gh/arviz-devs/arviz/pull/1753/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=arviz-devs#diff-YXJ2aXovZGF0YS9pb19lbWNlZS5weQ==) | `93.38% <0.00%> (-4.96%)` | :arrow_down: |
| [arviz/utils.py](https://codecov.io/gh/arviz-devs/arviz/pull/1753/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=arviz-devs#diff-YXJ2aXovdXRpbHMucHk=) | `88.81% <0.00%> (-0.33%)` | :arrow_down: |
| ... and [2 more](https://codecov.io/gh/arviz-devs/arviz/pull/1753/diff?src=pr&el=tree-more&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=arviz-devs) | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/arviz-devs/arviz/pull/1753?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=arviz-devs).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=arviz-devs)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/arviz-devs/arviz/pull/1753?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=arviz-devs). Last update [928f8f4...de5fdd1](https://codecov.io/gh/arviz-devs/arviz/pull/1753?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=arviz-devs). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=arviz-devs).

Rishabh261998: @sethaxen I have made the suggested changes, just the performance of the plot is the main issue now. Currently analysing each part and seeing if something can be optimized
sethaxen: @Rishabh261998 can you resolve the comments above that have been addressed so we know what is still open for discussion when we review the code?
Rishabh261998: @sethaxen have made all the changes. Please go through it. Also will add tests in some time.
michaelosthege: Any news on this PR? It would be great to see this in ArviZ!
A few months ago I tried to implement this method myself, but failed miserably.
sethaxen: Good question. @Rishabh261998, I'm sure you're busy with other things, but any chance you have the time to push this PR over the finish line? Else one of us can.
Rishabh261998: Hi @sethaxen, I am really sorry for not completing it.  As of now, I am a bit busy as I have joined a company and also in parallel, I am preparing for an exam. I think you guys can take it as I am not sure when I will get time for this. Thanks
sethaxen: > Hi @sethaxen, I am really sorry for not completing it. As of now, I am a bit busy as I have joined a company and also in parallel, I am preparing for an exam. I think you guys can take it as I am not sure when I will get time for this. Thanks

No problem at all, @Rishabh261998. Thanks to your excellent work on this, I don't think it will take us much to finish it up. Best of luck on the exam!

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 74bf414..702015c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,6 +9,7 @@
 * Add dict option to `from_cmdstan` log_likelihood parameter (as in `from_pystan`)
 * Unify model comparison API. Both `plot_compare`, `plot_elpd` can now take dicts of InferenceData or ELPDData ([1690](https://github.com/arviz-devs/arviz/pull/1690))
 * Change default for rcParam `stats.ic_pointwise` to True ([1690](https://github.com/arviz-devs/arviz/pull/1690))
+* Add new plot type: plot_ecdf ([1753](https://github.com/arviz-devs/arviz/pull/1753))
 
 
 ### Maintenance and fixes
diff --git a/arviz/plots/__init__.py b/arviz/plots/__init__.py
index f3c7ce0..8803c30 100644
--- a/arviz/plots/__init__.py
+++ b/arviz/plots/__init__.py
@@ -6,6 +6,7 @@ from .densityplot import plot_density
 from .distcomparisonplot import plot_dist_comparison
 from .distplot import plot_dist
 from .dotplot import plot_dot
+from .ecdfplot import plot_ecdf
 from .elpdplot import plot_elpd
 from .energyplot import plot_energy
 from .essplot import plot_ess
@@ -33,6 +34,7 @@ __all__ = [
     "plot_density",
     "plot_dist",
     "plot_dot",
+    "plot_ecdf",
     "plot_elpd",
     "plot_energy",
     "plot_ess",
diff --git a/arviz/plots/backends/bokeh/ecdfplot.py b/arviz/plots/backends/bokeh/ecdfplot.py
new file mode 100644
index 0000000..d756f62
--- /dev/null
+++ b/arviz/plots/backends/bokeh/ecdfplot.py
@@ -0,0 +1,73 @@
+"""Bokeh ecdfplot."""
+from matplotlib.colors import to_hex
+
+from ...plot_utils import _scale_fig_size
+from .. import show_layout
+from . import backend_kwarg_defaults, create_axes_grid
+
+
+def plot_ecdf(
+    x_coord,
+    y_coord,
+    x_bands,
+    lower,
+    higher,
+    confidence_bands,
+    plot_kwargs,
+    fill_kwargs,
+    plot_outline_kwargs,
+    figsize,
+    fill_band,
+    ax,
+    show,
+    backend_kwargs,
+):
+    """Bokeh ecdfplot."""
+    if backend_kwargs is None:
+        backend_kwargs = {}
+
+    backend_kwargs = {
+        **backend_kwarg_defaults(),
+        **backend_kwargs,
+    }
+
+    (figsize, *_) = _scale_fig_size(figsize, None)
+
+    if ax is None:
+        ax = create_axes_grid(
+            1,
+            figsize=figsize,
+            squeeze=True,
+            backend_kwargs=backend_kwargs,
+        )
+
+    if plot_kwargs is None:
+        plot_kwargs = {}
+
+    plot_kwargs.setdefault("mode", "after")
+
+    if fill_band:
+        if fill_kwargs is None:
+            fill_kwargs = {}
+        fill_kwargs.setdefault("fill_color", to_hex("C0"))
+        fill_kwargs.setdefault("fill_alpha", 0.2)
+    else:
+        if plot_outline_kwargs is None:
+            plot_outline_kwargs = {}
+        plot_outline_kwargs.setdefault("color", to_hex("C0"))
+        plot_outline_kwargs.setdefault("alpha", 0.2)
+
+    if confidence_bands:
+        ax.step(x_coord, y_coord, **plot_kwargs)
+
+        if fill_band:
+            ax.varea(x_bands, lower, higher, **fill_kwargs)
+        else:
+            ax.line(x_bands, lower, **plot_outline_kwargs)
+            ax.line(x_bands, higher, **plot_outline_kwargs)
+    else:
+        ax.step(x_coord, y_coord, **plot_kwargs)
+
+    show_layout(ax, show)
+
+    return ax
diff --git a/arviz/plots/backends/matplotlib/ecdfplot.py b/arviz/plots/backends/matplotlib/ecdfplot.py
new file mode 100644
index 0000000..8b992d4
--- /dev/null
+++ b/arviz/plots/backends/matplotlib/ecdfplot.py
@@ -0,0 +1,70 @@
+"""Matplotlib ecdfplot."""
+import matplotlib.pyplot as plt
+from matplotlib.colors import to_hex
+
+from ...plot_utils import _scale_fig_size
+from . import backend_kwarg_defaults, create_axes_grid, backend_show
+
+
+def plot_ecdf(
+    x_coord,
+    y_coord,
+    x_bands,
+    lower,
+    higher,
+    confidence_bands,
+    plot_kwargs,
+    fill_kwargs,
+    plot_outline_kwargs,
+    figsize,
+    fill_band,
+    ax,
+    show,
+    backend_kwargs,
+):
+    """Matplotlib ecdfplot."""
+    if backend_kwargs is None:
+        backend_kwargs = {}
+
+    backend_kwargs = {
+        **backend_kwarg_defaults(),
+        **backend_kwargs,
+    }
+
+    (figsize, _, _, _, _, _) = _scale_fig_size(figsize, None)
+    backend_kwargs.setdefault("figsize", figsize)
+    backend_kwargs["squeeze"] = True
+
+    if ax is None:
+        _, ax = create_axes_grid(1, backend_kwargs=backend_kwargs)
+
+    if plot_kwargs is None:
+        plot_kwargs = {}
+
+    plot_kwargs.setdefault("where", "post")
+
+    if fill_band:
+        if fill_kwargs is None:
+            fill_kwargs = {}
+        fill_kwargs.setdefault("step", "post")
+        fill_kwargs.setdefault("color", to_hex("C0"))
+        fill_kwargs.setdefault("alpha", 0.2)
+    else:
+        if plot_outline_kwargs is None:
+            plot_outline_kwargs = {}
+        plot_outline_kwargs.setdefault("where", "post")
+        plot_outline_kwargs.setdefault("color", to_hex("C0"))
+        plot_outline_kwargs.setdefault("alpha", 0.2)
+
+    ax.step(x_coord, y_coord, **plot_kwargs)
+
+    if confidence_bands:
+        if fill_band:
+            ax.fill_between(x_bands, lower, higher, **fill_kwargs)
+        else:
+            ax.plot(x_bands, lower, x_bands, higher, **plot_outline_kwargs)
+
+    if backend_show(show):
+        plt.show()
+
+    return ax
diff --git a/arviz/plots/ecdfplot.py b/arviz/plots/ecdfplot.py
new file mode 100644
index 0000000..5da1101
--- /dev/null
+++ b/arviz/plots/ecdfplot.py
@@ -0,0 +1,318 @@
+"""Plot ecdf or ecdf-difference plot with confidence bands."""
+import numpy as np
+from scipy.stats import uniform, binom
+
+from ..rcparams import rcParams
+from .plot_utils import get_plotting_function
+
+
+def plot_ecdf(
+    values,
+    values2=None,
+    cdf=None,
+    difference=False,
+    pit=False,
+    confidence_bands=None,
+    pointwise=False,
+    npoints=100,
+    num_trials=500,
+    fpr=0.05,
+    figsize=None,
+    fill_band=True,
+    plot_kwargs=None,
+    fill_kwargs=None,
+    plot_outline_kwargs=None,
+    ax=None,
+    show=None,
+    backend=None,
+    backend_kwargs=None,
+    **kwargs
+):
+    """Plot ECDF or ECDF-Difference Plot with Confidence bands.
+
+    This plot uses the simulated based algorithm presented in the paper "Graphical Test for
+    Discrete Uniformity and its Applications in Goodness of Fit Evaluation and
+    Multiple Sample Comparison" [1]_.
+
+    Parameters
+    ----------
+    values : array-like
+        Values to plot from an unknown continuous or discrete distribution
+    values2 : array-like, optional
+        Values to compare to the original sample
+    cdf : function, optional
+        Cumulative distribution function of the distribution to compare the original sample to
+    difference : bool, optional, Defaults False
+        If true then plot ECDF-difference plot otherwise ECDF plot
+    pit : bool, optional
+        If True plots the ECDF or ECDF-diff of PIT of sample
+    confidence_bands : bool, optional, Defaults True
+        If True plots the simultaneous or pointwise confidence bands with 1 - fpr confidence level
+    pointwise : bool, optional, Defaults False
+        If True plots pointwise confidence bands otherwise simultaneous bands
+    npoints : int, optional, Defaults 100
+        This denotes the granularity size of our plot
+        i.e the number of evaluation points for our ecdf or ecdf-difference plot
+    num_trials : int, optional, Defaults 500
+        The number of random ECDFs to generate to construct simultaneous confidence bands
+    fpr : float, optional, Defaults 0.05
+        The type I error rate s.t 1 - fpr denotes the confidence level of bands
+    figsize : tuple, optional
+        Figure size. If None it will be defined automatically.
+    fill_band : bool, optional
+        Use fill_between to mark the area inside the credible interval.
+        Otherwise, plot the border lines.
+    plot_kwargs : dict, optional
+        Additional kwargs passed to :func:`mpl:matplotlib.pyplot.step` or
+        :meth:`bokeh:bokeh.plotting.Figure.step`
+    fill_kwargs : dict, optional
+        Additional kwargs passed to :func:`mpl:matplotlib.pyplot.fill_between` or
+        :meth:`bokeh:bokeh.plotting.Figure.varea`
+    plot_outline_kwargs : dict, optional
+        Additional kwargs passed to :meth:`mpl:matplotlib.axes.Axes.plot` or
+        :meth:`bokeh:bokeh.plotting.Figure.line`
+    ax : axes, optional
+        Matplotlib axes or bokeh figures.
+    show : bool, optional
+        Call backend show function.
+    backend : str, optional
+        Select plotting backend {"matplotlib","bokeh"}. Default "matplotlib".
+    backend_kwargs : dict, optional
+        These are kwargs specific to the backend being used, passed to
+        :func:`mpl:matplotlib.pyplot.subplots` or
+        :meth:`bokeh:bokeh.plotting.figure`.
+
+    Returns
+    -------
+    axes : matplotlib axes or bokeh figures
+
+    References
+    ----------
+    .. [1] Säilynoja, T., Bürkner, P.C. and Vehtari, A., 2021. Graphical Test for
+        Discrete Uniformity and its Applications in Goodness of Fit Evaluation and
+        Multiple Sample Comparison. arXiv preprint arXiv:2103.10522.
+
+    Examples
+    --------
+    Plot ecdf plot for a given sample
+
+    .. plot::
+        :context: close-figs
+
+        >>> import arviz as az
+        >>> from scipy.stats import uniform, binom, norm
+
+        >>> sample = norm(0,1).rvs(1000)
+        >>> az.plot_ecdf(sample)
+
+    Plot ecdf plot with confidence bands for comparing a given sample w.r.t a given distribution
+
+    .. plot::
+        :context: close-figs
+
+        >>> distribution = norm(0,1)
+        >>> az.plot_ecdf(sample, cdf = distribution.cdf, confidence_bands = True)
+
+    Plot ecdf-difference plot with confidence bands for comparing a given sample
+    w.r.t a given distribution
+
+    .. plot::
+        :context: close-figs
+
+        >>> az.plot_ecdf(sample, cdf = distribution.cdf,
+            confidence_bands = True, difference = True)
+
+    Plot ecdf plot with confidence bands for PIT of sample for comparing a given sample
+    w.r.t a given distribution
+
+    .. plot::
+        :context: close-figs
+
+        >>> az.plot_ecdf(sample, cdf = distribution.cdf,
+            confidence_bands = True, pit = True)
+
+    Plot ecdf-difference plot with confidence bands for PIT of sample for comparing a given
+    sample w.r.t a given distribution
+
+    .. plot::
+        :context: close-figs
+
+        >>> az.plot_ecdf(sample, cdf = distribution.cdf,
+            confidence_bands = True, difference = True, pit = True)
+
+    You could also plot the above w.r.t another sample rather than a given distribution.
+    For eg: Plot ecdf-difference plot with confidence bands for PIT of sample for
+    comparing a given sample w.r.t a given sample
+
+    .. plot::
+        :context: close-figs
+
+        >>> sample2 = norm(0,1).rvs(5000)
+        >>> az.plot_ecdf(sample, sample2, confidence_bands = True, difference = True, pit = True)
+
+    """
+    if confidence_bands is None:
+        confidence_bands = (values2 is not None) or (cdf is not None)
+
+    if values2 is None and cdf is None and confidence_bands is True:
+        raise ValueError("For confidence bands you need to specify values2 or the cdf")
+
+    if cdf is not None and values2 is not None:
+        raise ValueError("To compare sample you need either cdf or values2 and not both")
+
+    if values2 is None and cdf is None and pit is True:
+        raise ValueError("For PIT specify either cdf or values2")
+
+    if values2 is None and cdf is None and difference is True:
+        raise ValueError("For ECDF difference plot need either cdf or values2")
+
+    if values2 is not None:
+        values2 = np.ravel(values2)
+        values2.sort()
+
+    values = np.ravel(values)
+    values.sort()
+
+    n = len(values)  # number of samples
+    ## This block computes gamma and uses it to get the upper and lower confidence bands
+    ## Here we check if we want confidence bands or not
+    if confidence_bands:
+        ## If plotting PIT then we find the PIT values of sample.
+        ## Basically here we generate the evaluation points(x) and find the PIT values.
+        ## z is the evaluation point for our uniform distribution in compute_gamma()
+        if pit:
+            x = np.linspace(1 / npoints, 1, npoints)
+            z = x
+            ## Finding PIT for our sample
+            probs = cdf(values) if cdf else compute_ecdf(values2, values) / len(values2)
+        else:
+            ## If not PIT use sample for plots and for evaluation points(x) use equally spaced
+            ## points between minimum and maximum of sample
+            ## For z we have used cdf(x)
+            x = np.linspace(values[0], values[-1], npoints)
+            z = cdf(x) if cdf else compute_ecdf(values2, x)
+            probs = values
+
+        ## Computing gamma
+        if not pointwise:
+            gamma = compute_gamma(n, z, npoints, num_trials, fpr)
+        else:
+            gamma = fpr
+        ## Using gamma to get the confidence intervals
+        lower, higher = get_lims(gamma, n, z)
+
+        ## This block is for whether to plot ECDF or ECDF-difference
+        if not difference:
+            ## We store the coordinates of our ecdf in x_coord, y_coord
+            x_coord, y_coord = get_ecdf_points(x, probs, difference)
+        else:
+            ## Here we subtract the ecdf value as here we are plotting the ECDF-difference
+            x_coord, y_coord = get_ecdf_points(x, probs, difference)
+            for i, x_i in enumerate(x):
+                y_coord[i] = y_coord[i] - (
+                    x_i if pit else cdf(x_i) if cdf else compute_ecdf(values2, x_i)
+                )
+
+            ## Similarly we subtract from the upper and lower bounds
+            if pit:
+                lower = lower - x
+                higher = higher - x
+            else:
+                lower = lower - (cdf(x) if cdf else compute_ecdf(values2, x))
+                higher = higher - (cdf(x) if cdf else compute_ecdf(values2, x))
+
+    else:
+        if pit:
+            x = np.linspace(1 / npoints, 1, npoints)
+            probs = cdf(values)
+        else:
+            x = np.linspace(values[0], values[-1], npoints)
+            probs = values
+
+        lower, higher = None, None
+        ## This block is for whether to plot ECDF or ECDF-difference
+        if not difference:
+            x_coord, y_coord = get_ecdf_points(x, probs, difference)
+        else:
+            ## Here we subtract the ecdf value as here we are plotting the ECDF-difference
+            x_coord, y_coord = get_ecdf_points(x, probs, difference)
+            for i, x_i in enumerate(x):
+                y_coord[i] = y_coord[i] - (
+                    x_i if pit else cdf(x_i) if cdf else compute_ecdf(values2, x_i)
+                )
+
+    ecdf_plot_args = dict(
+        x_coord=x_coord,
+        y_coord=y_coord,
+        x_bands=x,
+        lower=lower,
+        higher=higher,
+        confidence_bands=confidence_bands,
+        figsize=figsize,
+        fill_band=fill_band,
+        plot_kwargs=plot_kwargs,
+        fill_kwargs=fill_kwargs,
+        plot_outline_kwargs=plot_outline_kwargs,
+        ax=ax,
+        show=show,
+        backend_kwargs=backend_kwargs,
+        **kwargs
+    )
+
+    if backend is None:
+        backend = rcParams["plot.backend"]
+    backend = backend.lower()
+
+    plot = get_plotting_function("plot_ecdf", "ecdfplot", backend)
+    ax = plot(**ecdf_plot_args)
+
+    return ax
+
+
+def compute_ecdf(sample, z):
+    """Compute ECDF.
+
+    This function computes the ecdf value at the evaluation point
+        or a sorted set of evaluation points.
+    """
+    return np.searchsorted(sample, z, side="right") / len(sample)
+
+
+def get_ecdf_points(x, probs, difference):
+    """Compute the coordinates for the ecdf points using compute_ecdf."""
+    y = compute_ecdf(probs, x)
+
+    if not difference:
+        x = np.insert(x, 0, x[0])
+        y = np.insert(y, 0, 0)
+    return x, y
+
+
+def compute_gamma(n, z, npoints=None, num_trials=1000, fpr=0.05):
+    """Compute gamma for confidence interval calculation.
+
+    This function simulates an adjusted value of gamma to account for multiplicity
+    when forming an 1-fpr level confidence envelope for the ECDF of a sample.
+    """
+    if npoints is None:
+        npoints = n
+    gamma = []
+    for _ in range(num_trials):
+        unif_samples = uniform.rvs(0, 1, n)
+        unif_samples = np.sort(unif_samples)
+        gamma_m = 1000
+        ## Can compute ecdf for all the z together or one at a time.
+        f_z = compute_ecdf(unif_samples, z)
+        f_z = compute_ecdf(unif_samples, z)
+        gamma_m = 2 * min(
+            np.amin(binom.cdf(n * f_z, n, z)), np.amin(1 - binom.cdf(n * f_z - 1, n, z))
+        )
+        gamma.append(gamma_m)
+    return np.quantile(gamma, fpr)
+
+
+def get_lims(gamma, n, z):
+    """Compute the simultaneous 1 - fpr level confidence bands."""
+    lower = binom.ppf(gamma / 2, n, z)
+    upper = binom.ppf(1 - gamma / 2, n, z)
+    return lower / n, upper / n


Test patch:
diff --git a/arviz/tests/base_tests/test_plots_bokeh.py b/arviz/tests/base_tests/test_plots_bokeh.py
index 3fa7972..f9498ef 100644
--- a/arviz/tests/base_tests/test_plots_bokeh.py
+++ b/arviz/tests/base_tests/test_plots_bokeh.py
@@ -5,6 +5,7 @@ from copy import deepcopy
 import numpy as np
 import pytest
 from pandas import DataFrame  # pylint: disable=wrong-import-position
+from scipy.stats import norm  # pylint: disable=wrong-import-position
 
 from ...data import from_dict, load_arviz_data  # pylint: disable=wrong-import-position
 from ...plots import (  # pylint: disable=wrong-import-position
@@ -15,6 +16,7 @@ from ...plots import (  # pylint: disable=wrong-import-position
     plot_dist,
     plot_dist_comparison,
     plot_dot,
+    plot_ecdf,
     plot_elpd,
     plot_energy,
     plot_ess,
@@ -345,6 +347,26 @@ def test_plot_compare_no_ic(models):
     assert "['loo', 'waic']" in str(err.value)
 
 
+def test_plot_ecdf_basic():
+    data = np.random.randn(4, 1000)
+    axes = plot_ecdf(data, backend="bokeh", show=False)
+    assert axes is not None
+
+
+def test_plot_ecdf_values2():
+    data = np.random.randn(4, 1000)
+    data2 = np.random.randn(4, 500)
+    axes = plot_ecdf(data, data2, backend="bokeh", show=False)
+    assert axes is not None
+
+
+def test_plot_ecdf_cdf():
+    data = np.random.randn(4, 1000)
+    cdf = norm(0, 1).cdf
+    axes = plot_ecdf(data, cdf=cdf, backend="bokeh", show=False)
+    assert axes is not None
+
+
 @pytest.mark.parametrize(
     "kwargs", [{}, {"ic": "loo"}, {"xlabels": True, "scale": "log"}, {"threshold": 2}]
 )
diff --git a/arviz/tests/base_tests/test_plots_matplotlib.py b/arviz/tests/base_tests/test_plots_matplotlib.py
index dd4c5b5..b8464ec 100644
--- a/arviz/tests/base_tests/test_plots_matplotlib.py
+++ b/arviz/tests/base_tests/test_plots_matplotlib.py
@@ -8,7 +8,7 @@ import numpy as np
 import pytest
 from matplotlib import animation
 from pandas import DataFrame
-from scipy.stats import gaussian_kde
+from scipy.stats import gaussian_kde, norm
 
 from ...data import from_dict, load_arviz_data
 from ...plots import (
@@ -19,6 +19,7 @@ from ...plots import (
     plot_dist,
     plot_dist_comparison,
     plot_dot,
+    plot_ecdf,
     plot_elpd,
     plot_energy,
     plot_ess,
@@ -1191,6 +1192,26 @@ def test_kde_cumulative(limits):
     np.testing.assert_almost_equal(round(density[-1], 3), 1)
 
 
+def test_plot_ecdf_basic():
+    data = np.random.randn(4, 1000)
+    axes = plot_ecdf(data)
+    assert axes is not None
+
+
+def test_plot_ecdf_values2():
+    data = np.random.randn(4, 1000)
+    data2 = np.random.randn(4, 1000)
+    axes = plot_ecdf(data, data2)
+    assert axes is not None
+
+
+def test_plot_ecdf_cdf():
+    data = np.random.randn(4, 1000)
+    cdf = norm(0, 1).cdf
+    axes = plot_ecdf(data, cdf=cdf)
+    assert axes is not None
+
+
 @pytest.mark.parametrize(
     "kwargs",
     [