Plots posterior distributions and sample traces.

        Basically a wrapperfor pm.traceplot() plus some niceties, based partly on code
        from: https://pymc-devs.github.io/pymc3/notebooks/GLM-model-selection.html.

        Args:
            varnames (list): List of variable names to plot. If None, all
                eligible variables are plotted.
            ranefs (bool): If True (default), shows trace plots for
                individual random effects.
            transformed (bool): If False (default), excludes internally
                transformed variables from plotting.
            combined (bool): If True, concatenates all chains into one before
                plotting. If False (default), plots separately lines for
                each chain (on the same axes).
            hist (bool): If True, plots a histogram for each fixed effect,
                in addition to the kde plot. To prevent visual clutter,
                histograms are never plotted for random effects.
            bins (int): If hist is True, the number of bins in the histogram.
                Ignored if hist is False.
            kind (str): Either 'trace' (default) or 'priors'. If 'priors',
                this just internally calls Model.plot()

def plot(self, varnames=None, ranefs=True, transformed=False,
             combined=False, hist=False, bins=20, kind='trace'):
        '''Plots posterior distributions and sample traces.

        Basically a wrapperfor pm.traceplot() plus some niceties, based partly on code
        from: https://pymc-devs.github.io/pymc3/notebooks/GLM-model-selection.html.

        Args:
            varnames (list): List of variable names to plot. If None, all
                eligible variables are plotted.
            ranefs (bool): If True (default), shows trace plots for
                individual random effects.
            transformed (bool): If False (default), excludes internally
                transformed variables from plotting.
            combined (bool): If True, concatenates all chains into one before
                plotting. If False (default), plots separately lines for
                each chain (on the same axes).
            hist (bool): If True, plots a histogram for each fixed effect,
                in addition to the kde plot. To prevent visual clutter,
                histograms are never plotted for random effects.
            bins (int): If hist is True, the number of bins in the histogram.
                Ignored if hist is False.
            kind (str): Either 'trace' (default) or 'priors'. If 'priors',
                this just internally calls Model.plot()
        '''

        def _plot_row(data, row, title, hist=True):
            # density plot
            axes[row, 0].set_title(title)
            if pma is not None:
                arr = np.atleast_2d(data.values.T).T
                pma.kdeplot_op(axes[row, 0], arr, bw = 4.5)
            else:
                data.plot(kind='kde', ax=axes[row, 0], legend=False)

            # trace plot
            axes[row, 1].set_title(title)
            data.plot(kind='line', ax=axes[row, 1], legend=False)
            # histogram
            if hist:
                if pma is not None:
                    pma.histplot_op(axes[row, 0], arr)
                else:
                    data.plot(kind='hist', ax=axes[row, 0], legend=False,
                              normed=True, bins=bins)

        if kind == 'priors':
            return self.model.plot(varnames)

        # count the total number of rows in the plot
        names = self._filter_names(varnames, ranefs, transformed)
        random = [re.sub(r'_offset$', '', x) in self.model.random_terms
                  for x in names]
        rows = sum([len(self.level_dict[p]) if not r else 1
                    for p, r in zip(names, random)])

        # make the plot!
        fig, axes = plt.subplots(rows, 2, figsize=(12, 2*rows))
        if rows == 1:
            axes = np.array([axes])  # For consistent 2D indexing

        _select_args = {'ranefs': ranefs, 'transformed': transformed}

        # Create list of chains (for combined, just one list w/ all chains)
        chains = list(range(self.n_chains))
        if combined:
            chains = [chains]

        for c in chains:

            row = 0

            for p in names:

                df = self[p].to_df(chains=c, **_select_args)
                # if p == 'floor|county_sd':

                # fixed effects
                if re.sub(r'_offset$', '', p) not in self.model.random_terms:
                    for lev in self.level_dict[p]:
                        lev_df = df[lev]
                        _plot_row(lev_df, row, lev, hist)
                        row += 1

                # random effects
                else:
                    _plot_row(df, row, p, hist=False)  # too much clutter
                    row += 1

        fig.tight_layout()

        # For bernoulli models, tell user which event is being modeled
        if self.model.family.name == 'bernoulli':
            event = next(i for i, x in enumerate(self.model.y.data.flatten())
                         if x > .99)
            warnings.warn('Modeling the probability that {}==\'{}\''.format(
                self.model.y.name,
                str(self.model.clean_data[self.model.y.name][event])))

        return axes