def find_confidence_interval(x, pdf, confidence_level):
    return pdf[pdf > x].sum() - confidence_level

def density_contour(xdata, ydata, binsx, binsy, ax=None, levels_confidence=[0.68, 0.95, 0.99], range=None, **contour_kwargs):
    """ Create a density contour plot.
    Parameters
    ----------
    xdata : numpy.ndarray
    ydata : numpy.ndarray
    binsx : int
        Number of bins along x dimension
    binsy : int
        Number of bins along y dimension
    ax : matplotlib.Axes (optional)
        If supplied, plot the contour to this axis. Otherwise, open a new figure
    contour_kwargs : dict
        kwargs to be passed to pyplot.contour()
    """
    import scipy.optimize as so
    # nbins_x = len(binsx) - 1
    # nbins_y = len(binsy) - 1
    H, xedges, yedges = np.histogram2d(xdata, ydata, bins=[binsx, binsy], range=range, normed=True)
    x_bin_sizes = (xedges[1:] - xedges[:-1])
    y_bin_sizes = (yedges[1:] - yedges[:-1])
    pdf = (H * (x_bin_sizes * y_bin_sizes))
    levels = [so.brentq(find_confidence_interval, 0., 1., args=(pdf, lvl)) for lvl in levels_confidence]
    # one_sigma = so.brentq(find_confidence_interval, 0., 1., args = (pdf, 0.68))
    # two_sigma = so.brentq(find_confidence_interval, 0., 1., args = (pdf, 0.95))
    # three_sigma = so.brentq(find_confidence_interval, 0., 1., args = (pdf, 0.99))
    # levels = [one_sigma, two_sigma, three_sigma]
    X, Y = 0.5 * (xedges[1:] + xedges[:-1]), 0.5 * (yedges[1:] + yedges[:-1])
    Z = pdf.T
    if ax is None:
        contour = plt.tricontour(X, Y, Z, levels=levels[::-1], origin="lower", **contour_kwargs)
    else:
        # contour = ax.contour(X, Y, Z, levels=levels, origin="lower", **contour_kwargs)
        contour = ax.contour(X, Y, Z, levels=levels[::-1], origin="lower", **contour_kwargs)
    return contour