We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
Passing axes to plot_density fails with several datasets
**Describe the bug**
I would like to pass existing axes to plot_density with the `ax` keyword. This works fine if I plot data from a single dataset. But I get an error when I try to plot several datasets at once (I think the error arises when it tries to produce a legend).
Plotting several datasets without providing axes also works fine.

**To Reproduce**
```python
import arviz
import matplotlib.pyplot as plt

test_data = arviz.load_arviz_data('centered_eight')

# This works
fig, ax = plt.subplots(3, 3)
ax1 = arviz.plot_density(data=[test_data.posterior], var_names=['mu', 'theta'], ax=ax);
# This works as well
ax2 = arviz.plot_density(data=[test_data.prior, test_data.posterior], var_names=['mu', 'theta']);
# This does not work
fig3, ax3 = plt.subplots(3, 3)
arviz.plot_density(data=[test_data.prior, test_data.posterior], var_names=['mu', 'theta'], ax=ax3);
```

This is the error I get:
```python
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-133-7feead268671> in <module>
     11 # This does not work
     12 fig3, ax3 = plt.subplots(3, 3)
---> 13 arviz.plot_density(data=[test_data.prior, test_data.posterior], var_names=['mu', 'theta'], ax=ax3);

~/miniconda3/envs/sunode/lib/python3.7/site-packages/arviz/plots/densityplot.py in plot_density(data, group, data_labels, var_names, transform, credible_interval, point_estimate, colors, outline, hpd_markers, shade, bw, figsize, textsize, ax, backend, backend_kwargs, show)
    263     # TODO: Add backend kwargs
    264     plot = get_plotting_function("plot_density", "densityplot", backend)
--> 265     ax = plot(**plot_density_kwargs)
    266     return ax

~/miniconda3/envs/sunode/lib/python3.7/site-packages/arviz/plots/backends/matplotlib/densityplot.py in plot_density(ax, all_labels, to_plot, colors, bw, figsize, length_plotters, rows, cols, titlesize, xt_labelsize, linewidth, markersize, credible_interval, point_estimate, hpd_markers, outline, shade, n_data, data_labels, backend_kwargs, show)
     75     if n_data > 1:
     76         for m_idx, label in enumerate(data_labels):
---> 77             ax[0].plot([], label=label, c=colors[m_idx], markersize=markersize)
     78         ax[0].legend(fontsize=xt_labelsize)
     79 

AttributeError: 'numpy.ndarray' object has no attribute 'plot'
```

**Additional context**
ArviZ version: 0.7.0
matplotlib version: 3.2.1


--- END ISSUE ---

Below are some code segments, each from a relevant file. One or more of these files may contain bugs.
--- BEGIN FILES ---
Path: `arviz/plots/backends/matplotlib/densityplot.py`
Content:
```
1 """Matplotlib Densityplot."""
2 import matplotlib.pyplot as plt
3 import numpy as np
4 
5 from . import backend_show
6 from ....stats import hdi
7 from ...plot_utils import (
8     make_label,
9     _create_axes_grid,
10     calculate_point_estimate,
11 )
12 from ....numeric_utils import _fast_kde, get_bins
13 
14 
15 def plot_density(
16     ax,
17     all_labels,
18     to_plot,
19     colors,
20     bw,
21     figsize,
22     length_plotters,
23     rows,
24     cols,
25     titlesize,
26     xt_labelsize,
27     linewidth,
28     markersize,
29     hdi_prob,
30     point_estimate,
31     hdi_markers,
32     outline,
33     shade,
34     n_data,
35     data_labels,
36     backend_kwargs,
37     show,
38 ):
39     """Matplotlib densityplot."""
40     if ax is None:
41         _, ax = _create_axes_grid(
42             length_plotters,
43             rows,
44             cols,
45             figsize=figsize,
46             squeeze=False,
47             backend="matplotlib",
48             backend_kwargs=backend_kwargs,
49         )
50     else:
51         ax = np.atleast_2d(ax)
52 
53     axis_map = {label: ax_ for label, ax_ in zip(all_labels, np.ravel(ax))}
54 
55     for m_idx, plotters in enumerate(to_plot):
56         for var_name, selection, values in plotters:
57             label = make_label(var_name, selection)
58             _d_helper(
59                 values.flatten(),
60                 label,
61                 colors[m_idx],
62                 bw,
63                 titlesize,
64                 xt_labelsize,
65                 linewidth,
66                 markersize,
67                 hdi_prob,
68                 point_estimate,
69                 hdi_markers,
70                 outline,
71                 shade,
72                 axis_map[label],
73             )
74 
75     if n_data > 1:
76         for m_idx, label in enumerate(data_labels):
77             ax[0].plot([], label=label, c=colors[m_idx], markersize=markersize)
78         ax[0].legend(fontsize=xt_labelsize)
79 
80     if backend_show(show):
81         plt.show()
82 
83     return ax
84 
85 
86 def _d_helper(
87     vec,
88     vname,
89     color,
90     bw,
91     titlesize,
92     xt_labelsize,
93     linewidth,
94     markersize,
95     hdi_prob,
96     point_estimate,
97     hdi_markers,
98     outline,
99     shade,
100     ax,
101 ):
102     """Plot an individual dimension.
103 
104     Parameters
105     ----------
106     vec : array
107         1D array from trace
108     vname : str
109         variable name
110     color : str
111         matplotlib color
112     bw : float
113         Bandwidth scaling factor. Should be larger than 0. The higher this number the smoother the
114         KDE will be. Defaults to 4.5 which is essentially the same as the Scott's rule of thumb
115         (the default used rule by SciPy).
116     titlesize : float
117         font size for title
118     xt_labelsize : float
119        fontsize for xticks
120     linewidth : float
121         Thickness of lines
122     markersize : float
123         Size of markers
124     hdi_prob : float
125         Probability for the highest density interval. Defaults to 0.94
126     point_estimate : Optional[str]
127         Plot point estimate per variable. Values should be 'mean', 'median', 'mode' or None.
128         Defaults to 'auto' i.e. it falls back to default set in rcParams.
129     shade : float
130         Alpha blending value for the shaded area under the curve, between 0 (no shade) and 1
131         (opaque). Defaults to 0.
132     ax : matplotlib axes
133     """
134     if vec.dtype.kind == "f":
135         if hdi_prob != 1:
136             hdi_ = hdi(vec, hdi_prob, multimodal=False)
137             new_vec = vec[(vec >= hdi_[0]) & (vec <= hdi_[1])]
138         else:
139             new_vec = vec
140 
141         density, xmin, xmax = _fast_kde(new_vec, bw=bw)
142         density *= hdi_prob
143         x = np.linspace(xmin, xmax, len(density))
144         ymin = density[0]
145         ymax = density[-1]
146 
147         if outline:
148             ax.plot(x, density, color=color, lw=linewidth)
149             ax.plot([xmin, xmin], [-ymin / 100, ymin], color=color, ls="-", lw=linewidth)
150             ax.plot([xmax, xmax], [-ymax / 100, ymax], color=color, ls="-", lw=linewidth)
151 
152         if shade:
153             ax.fill_between(x, density, color=color, alpha=shade)
154 
155     else:
156         xmin, xmax = hdi(vec, hdi_prob, multimodal=False)
157         bins = get_bins(vec)
158         if outline:
159             ax.hist(vec, bins=bins, color=color, histtype="step", align="left")
160         if shade:
161             ax.hist(vec, bins=bins, color=color, alpha=shade)
162 
163     if hdi_markers:
164         ax.plot(xmin, 0, hdi_markers, color=color, markeredgecolor="k", markersize=markersize)
165         ax.plot(xmax, 0, hdi_markers, color=color, markeredgecolor="k", markersize=markersize)
166 
167     if point_estimate is not None:
168         est = calculate_point_estimate(point_estimate, vec, bw)
169         ax.plot(est, 0, "o", color=color, markeredgecolor="k", markersize=markersize)
170 
171     ax.set_yticks([])
172     ax.set_title(vname, fontsize=titlesize, wrap=True)
173     for pos in ["left", "right", "top"]:
174         ax.spines[pos].set_visible(False)
175     ax.tick_params(labelsize=xt_labelsize)
176
```
--- END FILES ---

Please first localize the bug based on the issue statement, and then generate a patch according to the `git diff` format fenced by three backticks.

Here is an example:

```diff
diff --git a/examples/server_async.py b/examples/server_async.py
--- a/examples/server_async.py
+++ b/examples/server_async.py
@@ -313,4 +313,4 @@
 
 
 if __name__ == "__main__":
-    asyncio.run(run_async_server("."), debug=True)
+    asyncio.run(run_async_server(), debug=True)
diff --git a/examples/server_sync.py b/examples/server_sync.py
--- a/examples/server_sync.py
+++ b/examples/server_sync.py
@@ -313,5 +313,5 @@
 
 
 if __name__ == "__main__":
-    server = run_sync_server(".")
+    server = run_sync_server()
     server.shutdown()

```

diff --git a/arviz/plots/backends/matplotlib/densityplot.py b/arviz/plots/backends/matplotlib/densityplot.py
--- a/arviz/plots/backends/matplotlib/densityplot.py
+++ b/arviz/plots/backends/matplotlib/densityplot.py
@@ -74,8 +74,8 @@
 
     if n_data > 1:
         for m_idx, label in enumerate(data_labels):
-            ax[0].plot([], label=label, c=colors[m_idx], markersize=markersize)
-        ax[0].legend(fontsize=xt_labelsize)
+            ax.item(0).plot([], label=label, c=colors[m_idx], markersize=markersize)
+        ax.item(0).legend(fontsize=xt_labelsize)
 
     if backend_show(show):
         plt.show()

{"golden_diff": "diff --git a/arviz/plots/backends/matplotlib/densityplot.py b/arviz/plots/backends/matplotlib/densityplot.py\n--- a/arviz/plots/backends/matplotlib/densityplot.py\n+++ b/arviz/plots/backends/matplotlib/densityplot.py\n@@ -74,8 +74,8 @@\n \n     if n_data > 1:\n         for m_idx, label in enumerate(data_labels):\n-            ax[0].plot([], label=label, c=colors[m_idx], markersize=markersize)\n-        ax[0].legend(fontsize=xt_labelsize)\n+            ax.item(0).plot([], label=label, c=colors[m_idx], markersize=markersize)\n+        ax.item(0).legend(fontsize=xt_labelsize)\n \n     if backend_show(show):\n         plt.show()\n", "issue": "Passing axes to plot_density fails with several datasets\n**Describe the bug**\r\nI would like to pass existing axes to plot_density with the `ax` keyword. This works fine if I plot data from a single dataset. But I get an error when I try to plot several datasets at once (I think the error arises when it tries to produce a legend).\r\nPlotting several datasets without providing axes also works fine.\r\n\r\n**To Reproduce**\r\n```python\r\nimport arviz\r\nimport matplotlib.pyplot as plt\r\n\r\ntest_data = arviz.load_arviz_data('centered_eight')\r\n\r\n# This works\r\nfig, ax = plt.subplots(3, 3)\r\nax1 = arviz.plot_density(data=[test_data.posterior], var_names=['mu', 'theta'], ax=ax);\r\n# This works as well\r\nax2 = arviz.plot_density(data=[test_data.prior, test_data.posterior], var_names=['mu', 'theta']);\r\n# This does not work\r\nfig3, ax3 = plt.subplots(3, 3)\r\narviz.plot_density(data=[test_data.prior, test_data.posterior], var_names=['mu', 'theta'], ax=ax3);\r\n```\r\n\r\nThis is the error I get:\r\n```python\r\n---------------------------------------------------------------------------\r\nAttributeError                            Traceback (most recent call last)\r\n<ipython-input-133-7feead268671> in <module>\r\n     11 # This does not work\r\n     12 fig3, ax3 = plt.subplots(3, 3)\r\n---> 13 arviz.plot_density(data=[test_data.prior, test_data.posterior], var_names=['mu', 'theta'], ax=ax3);\r\n\r\n~/miniconda3/envs/sunode/lib/python3.7/site-packages/arviz/plots/densityplot.py in plot_density(data, group, data_labels, var_names, transform, credible_interval, point_estimate, colors, outline, hpd_markers, shade, bw, figsize, textsize, ax, backend, backend_kwargs, show)\r\n    263     # TODO: Add backend kwargs\r\n    264     plot = get_plotting_function(\"plot_density\", \"densityplot\", backend)\r\n--> 265     ax = plot(**plot_density_kwargs)\r\n    266     return ax\r\n\r\n~/miniconda3/envs/sunode/lib/python3.7/site-packages/arviz/plots/backends/matplotlib/densityplot.py in plot_density(ax, all_labels, to_plot, colors, bw, figsize, length_plotters, rows, cols, titlesize, xt_labelsize, linewidth, markersize, credible_interval, point_estimate, hpd_markers, outline, shade, n_data, data_labels, backend_kwargs, show)\r\n     75     if n_data > 1:\r\n     76         for m_idx, label in enumerate(data_labels):\r\n---> 77             ax[0].plot([], label=label, c=colors[m_idx], markersize=markersize)\r\n     78         ax[0].legend(fontsize=xt_labelsize)\r\n     79 \r\n\r\nAttributeError: 'numpy.ndarray' object has no attribute 'plot'\r\n```\r\n\r\n**Additional context**\r\nArviZ version: 0.7.0\r\nmatplotlib version: 3.2.1\r\n\n", "before_files": [{"content": "\"\"\"Matplotlib Densityplot.\"\"\"\nimport matplotlib.pyplot as plt\nimport numpy as np\n\nfrom . import backend_show\nfrom ....stats import hdi\nfrom ...plot_utils import (\n    make_label,\n    _create_axes_grid,\n    calculate_point_estimate,\n)\nfrom ....numeric_utils import _fast_kde, get_bins\n\n\ndef plot_density(\n    ax,\n    all_labels,\n    to_plot,\n    colors,\n    bw,\n    figsize,\n    length_plotters,\n    rows,\n    cols,\n    titlesize,\n    xt_labelsize,\n    linewidth,\n    markersize,\n    hdi_prob,\n    point_estimate,\n    hdi_markers,\n    outline,\n    shade,\n    n_data,\n    data_labels,\n    backend_kwargs,\n    show,\n):\n    \"\"\"Matplotlib densityplot.\"\"\"\n    if ax is None:\n        _, ax = _create_axes_grid(\n            length_plotters,\n            rows,\n            cols,\n            figsize=figsize,\n            squeeze=False,\n            backend=\"matplotlib\",\n            backend_kwargs=backend_kwargs,\n        )\n    else:\n        ax = np.atleast_2d(ax)\n\n    axis_map = {label: ax_ for label, ax_ in zip(all_labels, np.ravel(ax))}\n\n    for m_idx, plotters in enumerate(to_plot):\n        for var_name, selection, values in plotters:\n            label = make_label(var_name, selection)\n            _d_helper(\n                values.flatten(),\n                label,\n                colors[m_idx],\n                bw,\n                titlesize,\n                xt_labelsize,\n                linewidth,\n                markersize,\n                hdi_prob,\n                point_estimate,\n                hdi_markers,\n                outline,\n                shade,\n                axis_map[label],\n            )\n\n    if n_data > 1:\n        for m_idx, label in enumerate(data_labels):\n            ax[0].plot([], label=label, c=colors[m_idx], markersize=markersize)\n        ax[0].legend(fontsize=xt_labelsize)\n\n    if backend_show(show):\n        plt.show()\n\n    return ax\n\n\ndef _d_helper(\n    vec,\n    vname,\n    color,\n    bw,\n    titlesize,\n    xt_labelsize,\n    linewidth,\n    markersize,\n    hdi_prob,\n    point_estimate,\n    hdi_markers,\n    outline,\n    shade,\n    ax,\n):\n    \"\"\"Plot an individual dimension.\n\n    Parameters\n    ----------\n    vec : array\n        1D array from trace\n    vname : str\n        variable name\n    color : str\n        matplotlib color\n    bw : float\n        Bandwidth scaling factor. Should be larger than 0. The higher this number the smoother the\n        KDE will be. Defaults to 4.5 which is essentially the same as the Scott's rule of thumb\n        (the default used rule by SciPy).\n    titlesize : float\n        font size for title\n    xt_labelsize : float\n       fontsize for xticks\n    linewidth : float\n        Thickness of lines\n    markersize : float\n        Size of markers\n    hdi_prob : float\n        Probability for the highest density interval. Defaults to 0.94\n    point_estimate : Optional[str]\n        Plot point estimate per variable. Values should be 'mean', 'median', 'mode' or None.\n        Defaults to 'auto' i.e. it falls back to default set in rcParams.\n    shade : float\n        Alpha blending value for the shaded area under the curve, between 0 (no shade) and 1\n        (opaque). Defaults to 0.\n    ax : matplotlib axes\n    \"\"\"\n    if vec.dtype.kind == \"f\":\n        if hdi_prob != 1:\n            hdi_ = hdi(vec, hdi_prob, multimodal=False)\n            new_vec = vec[(vec >= hdi_[0]) & (vec <= hdi_[1])]\n        else:\n            new_vec = vec\n\n        density, xmin, xmax = _fast_kde(new_vec, bw=bw)\n        density *= hdi_prob\n        x = np.linspace(xmin, xmax, len(density))\n        ymin = density[0]\n        ymax = density[-1]\n\n        if outline:\n            ax.plot(x, density, color=color, lw=linewidth)\n            ax.plot([xmin, xmin], [-ymin / 100, ymin], color=color, ls=\"-\", lw=linewidth)\n            ax.plot([xmax, xmax], [-ymax / 100, ymax], color=color, ls=\"-\", lw=linewidth)\n\n        if shade:\n            ax.fill_between(x, density, color=color, alpha=shade)\n\n    else:\n        xmin, xmax = hdi(vec, hdi_prob, multimodal=False)\n        bins = get_bins(vec)\n        if outline:\n            ax.hist(vec, bins=bins, color=color, histtype=\"step\", align=\"left\")\n        if shade:\n            ax.hist(vec, bins=bins, color=color, alpha=shade)\n\n    if hdi_markers:\n        ax.plot(xmin, 0, hdi_markers, color=color, markeredgecolor=\"k\", markersize=markersize)\n        ax.plot(xmax, 0, hdi_markers, color=color, markeredgecolor=\"k\", markersize=markersize)\n\n    if point_estimate is not None:\n        est = calculate_point_estimate(point_estimate, vec, bw)\n        ax.plot(est, 0, \"o\", color=color, markeredgecolor=\"k\", markersize=markersize)\n\n    ax.set_yticks([])\n    ax.set_title(vname, fontsize=titlesize, wrap=True)\n    for pos in [\"left\", \"right\", \"top\"]:\n        ax.spines[pos].set_visible(False)\n    ax.tick_params(labelsize=xt_labelsize)\n", "path": "arviz/plots/backends/matplotlib/densityplot.py"}], "after_files": [{"content": "\"\"\"Matplotlib Densityplot.\"\"\"\nimport matplotlib.pyplot as plt\nimport numpy as np\n\nfrom . import backend_show\nfrom ....stats import hdi\nfrom ...plot_utils import (\n    make_label,\n    _create_axes_grid,\n    calculate_point_estimate,\n)\nfrom ....numeric_utils import _fast_kde, get_bins\n\n\ndef plot_density(\n    ax,\n    all_labels,\n    to_plot,\n    colors,\n    bw,\n    figsize,\n    length_plotters,\n    rows,\n    cols,\n    titlesize,\n    xt_labelsize,\n    linewidth,\n    markersize,\n    hdi_prob,\n    point_estimate,\n    hdi_markers,\n    outline,\n    shade,\n    n_data,\n    data_labels,\n    backend_kwargs,\n    show,\n):\n    \"\"\"Matplotlib densityplot.\"\"\"\n    if ax is None:\n        _, ax = _create_axes_grid(\n            length_plotters,\n            rows,\n            cols,\n            figsize=figsize,\n            squeeze=False,\n            backend=\"matplotlib\",\n            backend_kwargs=backend_kwargs,\n        )\n    else:\n        ax = np.atleast_2d(ax)\n\n    axis_map = {label: ax_ for label, ax_ in zip(all_labels, np.ravel(ax))}\n\n    for m_idx, plotters in enumerate(to_plot):\n        for var_name, selection, values in plotters:\n            label = make_label(var_name, selection)\n            _d_helper(\n                values.flatten(),\n                label,\n                colors[m_idx],\n                bw,\n                titlesize,\n                xt_labelsize,\n                linewidth,\n                markersize,\n                hdi_prob,\n                point_estimate,\n                hdi_markers,\n                outline,\n                shade,\n                axis_map[label],\n            )\n\n    if n_data > 1:\n        for m_idx, label in enumerate(data_labels):\n            ax.item(0).plot([], label=label, c=colors[m_idx], markersize=markersize)\n        ax.item(0).legend(fontsize=xt_labelsize)\n\n    if backend_show(show):\n        plt.show()\n\n    return ax\n\n\ndef _d_helper(\n    vec,\n    vname,\n    color,\n    bw,\n    titlesize,\n    xt_labelsize,\n    linewidth,\n    markersize,\n    hdi_prob,\n    point_estimate,\n    hdi_markers,\n    outline,\n    shade,\n    ax,\n):\n    \"\"\"Plot an individual dimension.\n\n    Parameters\n    ----------\n    vec : array\n        1D array from trace\n    vname : str\n        variable name\n    color : str\n        matplotlib color\n    bw : float\n        Bandwidth scaling factor. Should be larger than 0. The higher this number the smoother the\n        KDE will be. Defaults to 4.5 which is essentially the same as the Scott's rule of thumb\n        (the default used rule by SciPy).\n    titlesize : float\n        font size for title\n    xt_labelsize : float\n       fontsize for xticks\n    linewidth : float\n        Thickness of lines\n    markersize : float\n        Size of markers\n    hdi_prob : float\n        Probability for the highest density interval. Defaults to 0.94\n    point_estimate : Optional[str]\n        Plot point estimate per variable. Values should be 'mean', 'median', 'mode' or None.\n        Defaults to 'auto' i.e. it falls back to default set in rcParams.\n    shade : float\n        Alpha blending value for the shaded area under the curve, between 0 (no shade) and 1\n        (opaque). Defaults to 0.\n    ax : matplotlib axes\n    \"\"\"\n    if vec.dtype.kind == \"f\":\n        if hdi_prob != 1:\n            hdi_ = hdi(vec, hdi_prob, multimodal=False)\n            new_vec = vec[(vec >= hdi_[0]) & (vec <= hdi_[1])]\n        else:\n            new_vec = vec\n\n        density, xmin, xmax = _fast_kde(new_vec, bw=bw)\n        density *= hdi_prob\n        x = np.linspace(xmin, xmax, len(density))\n        ymin = density[0]\n        ymax = density[-1]\n\n        if outline:\n            ax.plot(x, density, color=color, lw=linewidth)\n            ax.plot([xmin, xmin], [-ymin / 100, ymin], color=color, ls=\"-\", lw=linewidth)\n            ax.plot([xmax, xmax], [-ymax / 100, ymax], color=color, ls=\"-\", lw=linewidth)\n\n        if shade:\n            ax.fill_between(x, density, color=color, alpha=shade)\n\n    else:\n        xmin, xmax = hdi(vec, hdi_prob, multimodal=False)\n        bins = get_bins(vec)\n        if outline:\n            ax.hist(vec, bins=bins, color=color, histtype=\"step\", align=\"left\")\n        if shade:\n            ax.hist(vec, bins=bins, color=color, alpha=shade)\n\n    if hdi_markers:\n        ax.plot(xmin, 0, hdi_markers, color=color, markeredgecolor=\"k\", markersize=markersize)\n        ax.plot(xmax, 0, hdi_markers, color=color, markeredgecolor=\"k\", markersize=markersize)\n\n    if point_estimate is not None:\n        est = calculate_point_estimate(point_estimate, vec, bw)\n        ax.plot(est, 0, \"o\", color=color, markeredgecolor=\"k\", markersize=markersize)\n\n    ax.set_yticks([])\n    ax.set_title(vname, fontsize=titlesize, wrap=True)\n    for pos in [\"left\", \"right\", \"top\"]:\n        ax.spines[pos].set_visible(False)\n    ax.tick_params(labelsize=xt_labelsize)\n", "path": "arviz/plots/backends/matplotlib/densityplot.py"}]}