def freq_mag(magnitudes, completeness, max_mag, binsize=0.2, **kwargs):
    """
    Plot a frequency-magnitude histogram and cumulative density plot.

    Currently this will compute a b-value, for a given completeness.
    B-value is computed by linear fitting to section of curve between
    completeness and max_mag.

    :type magnitudes: list
    :param magnitudes: list of float of magnitudes
    :type completeness: float
    :param completeness: Level to compute the b-value above
    :type max_mag: float
    :param max_mag: Maximum magnitude to try and fit a b-value to
    :type binsize: float
    :param binsize: Width of histogram bins, defaults to 0.2

    :returns: :class:`matplotlib.figure.Figure`

    .. Note::
        See :func:`eqcorrscan.utils.mag_calc.calc_b_value` for a least-squares
        method of estimating completeness and b-value. For estimating maximum
        curvature see :func:`eqcorrscan.utils.mag_calc.calc_max_curv`.

    .. rubric:: Example

    >>> from obspy.clients.fdsn import Client
    >>> from obspy import UTCDateTime
    >>> from eqcorrscan.utils.plotting import freq_mag
    >>> client = Client('IRIS')
    >>> t1 = UTCDateTime('2012-03-26T00:00:00')
    >>> t2 = t1 + (3 * 86400)
    >>> catalog = client.get_events(starttime=t1, endtime=t2, minmagnitude=3)
    >>> magnitudes = [event.preferred_magnitude().mag for event in catalog]
    >>> freq_mag(magnitudes, completeness=4, max_mag=7) # doctest: +SKIP

    .. plot::

        from obspy.clients.fdsn import Client
        from obspy import UTCDateTime
        from eqcorrscan.utils.plotting import freq_mag
        client = Client('IRIS')
        t1 = UTCDateTime('2012-03-26T00:00:00')
        t2 = t1 + (3 * 86400)
        catalog = client.get_events(starttime=t1, endtime=t2, minmagnitude=3)
        magnitudes = [event.preferred_magnitude().mag for event in catalog]
        freq_mag(magnitudes, completeness=4, max_mag=7)
    """
    import matplotlib.pyplot as plt
    # Ensure magnitudes are sorted
    magnitudes.sort()
    # Check that there are no nans or infs
    if np.isnan(magnitudes).any():
        warnings.warn('Found nan values, removing them')
        magnitudes = [mag for mag in magnitudes if not np.isnan(mag)]
    if np.isinf(magnitudes).any():
        warnings.warn('Found inf values, removing them')
        magnitudes = [mag for mag in magnitudes if not np.isinf(mag)]
    fig, ax1 = plt.subplots()
    # Set up the bins, the bin-size could be a variables
    bins = np.arange(int(min(magnitudes) - 1), int(max(magnitudes) + 1),
                     binsize)
    n, bins, patches = ax1.hist(magnitudes, bins, facecolor='Black',
                                alpha=0.5, label='Magnitudes')
    ax1.set_ylabel('Frequency')
    ax1.set_ylim([0, max(n) + 0.5 * max(n)])
    plt.xlabel('Magnitude')
    # Now make the cumulative density function
    counts = Counter(magnitudes)
    cdf = np.zeros(len(counts))
    mag_steps = np.zeros(len(counts))
    for i, magnitude in enumerate(sorted(counts.keys(), reverse=True)):
        mag_steps[i] = magnitude
        if i > 0:
            cdf[i] = cdf[i - 1] + counts[magnitude]
        else:
            cdf[i] = counts[magnitude]
    ax2 = ax1.twinx()
    # ax2.scatter(magnitudes, np.log10(cdf), c='k', marker='+', s=20, lw=2,
    ax2.scatter(mag_steps, np.log10(cdf), c='k', marker='+', s=20, lw=2,
                label='Magnitude cumulative density')
    # Now we want to calculate the b-value and plot the fit
    x = []
    y = []
    for i, magnitude in enumerate(mag_steps):
        if magnitude >= completeness <= max_mag:
            x.append(magnitude)
            y.append(cdf[i])
    fit = np.polyfit(x, np.log10(y), 1)
    fit_fn = np.poly1d(fit)
    ax2.plot(magnitudes, fit_fn(magnitudes), '--k',
             label='GR trend, b-value = ' + str(abs(fit[0]))[0:4] +
             '\n $M_C$ = ' + str(completeness))
    ax2.set_ylabel('$Log_{10}$ of cumulative density')
    plt.xlim([min(magnitudes) - 0.1, max(magnitudes) + 0.2])
    plt.ylim([min(np.log10(cdf)) - 0.5, max(np.log10(cdf)) + 1.0])
    plt.legend(loc=2)
    fig = _finalise_figure(fig=fig, **kwargs)  # pragma: no cover
    return fig
Plot a frequency-magnitude histogram and cumulative density plot.

    Currently this will compute a b-value, for a given completeness.
    B-value is computed by linear fitting to section of curve between
    completeness and max_mag.

    :type magnitudes: list
    :param magnitudes: list of float of magnitudes
    :type completeness: float
    :param completeness: Level to compute the b-value above
    :type max_mag: float
    :param max_mag: Maximum magnitude to try and fit a b-value to
    :type binsize: float
    :param binsize: Width of histogram bins, defaults to 0.2

    :returns: :class:`matplotlib.figure.Figure`

    .. Note::
        See :func:`eqcorrscan.utils.mag_calc.calc_b_value` for a least-squares
        method of estimating completeness and b-value. For estimating maximum
        curvature see :func:`eqcorrscan.utils.mag_calc.calc_max_curv`.

    .. rubric:: Example

    >>> from obspy.clients.fdsn import Client
    >>> from obspy import UTCDateTime
    >>> from eqcorrscan.utils.plotting import freq_mag
    >>> client = Client('IRIS')
    >>> t1 = UTCDateTime('2012-03-26T00:00:00')
    >>> t2 = t1 + (3 * 86400)
    >>> catalog = client.get_events(starttime=t1, endtime=t2, minmagnitude=3)
    >>> magnitudes = [event.preferred_magnitude().mag for event in catalog]
    >>> freq_mag(magnitudes, completeness=4, max_mag=7) # doctest: +SKIP

    .. plot::

        from obspy.clients.fdsn import Client
        from obspy import UTCDateTime
        from eqcorrscan.utils.plotting import freq_mag
        client = Client('IRIS')
        t1 = UTCDateTime('2012-03-26T00:00:00')
        t2 = t1 + (3 * 86400)
        catalog = client.get_events(starttime=t1, endtime=t2, minmagnitude=3)
        magnitudes = [event.preferred_magnitude().mag for event in catalog]
        freq_mag(magnitudes, completeness=4, max_mag=7)