def calc_b_value(magnitudes, completeness, max_mag=None, plotvar=True):
b_values = []<EOL>counts = Counter(magnitudes)<EOL>cdf = np.zeros(len(counts))<EOL>mag_steps = np.zeros(len(counts))<EOL>for i, magnitude in enumerate(sorted(counts.keys(), reverse=True)):<EOL><INDENT>mag_steps[i] = magnitude<EOL>if i > <NUM_LIT:0>:<EOL><INDENT>cdf[i] = cdf[i - <NUM_LIT:1>] + counts[magnitude]<EOL><DEDENT>else:<EOL><INDENT>cdf[i] = counts[magnitude]<EOL><DEDENT><DEDENT>if not max_mag:<EOL><INDENT>max_mag = max(magnitudes)<EOL><DEDENT>for m_c in completeness:<EOL><INDENT>if m_c >= max_mag or m_c >= max(magnitudes):<EOL><INDENT>warnings.warn('<STR_LIT>' %<EOL>str(m_c))<EOL>break<EOL><DEDENT>complete_mags = []<EOL>complete_freq = []<EOL>for i, mag in enumerate(mag_steps):<EOL><INDENT>if mag >= m_c <= max_mag:<EOL><INDENT>complete_mags.append(mag)<EOL>complete_freq.append(np.log10(cdf[i]))<EOL><DEDENT><DEDENT>if len(complete_mags) < <NUM_LIT:4>:<EOL><INDENT>warnings.warn('<STR_LIT>' + str(m_c) +<EOL>'<STR_LIT>')<EOL>break<EOL><DEDENT>fit = np.polyfit(complete_mags, complete_freq, <NUM_LIT:1>, full=True)<EOL>predicted_freqs = [fit[<NUM_LIT:0>][<NUM_LIT:1>] - abs(fit[<NUM_LIT:0>][<NUM_LIT:0>] * M)<EOL>for M in complete_mags]<EOL>r = <NUM_LIT:100> - ((np.sum([abs(complete_freq[i] - predicted_freqs[i])<EOL>for i in range(len(complete_freq))]) * <NUM_LIT:100>) /<EOL>np.sum(complete_freq))<EOL>b_values.append((m_c, abs(fit[<NUM_LIT:0>][<NUM_LIT:0>]), r, str(len(complete_mags))))<EOL><DEDENT>if plotvar:<EOL><INDENT>fig, ax1 = plt.subplots()<EOL>b_vals = ax1.scatter(list(zip(*b_values))[<NUM_LIT:0>], list(zip(*b_values))[<NUM_LIT:1>],<EOL>c='<STR_LIT:k>')<EOL>resid = ax1.scatter(list(zip(*b_values))[<NUM_LIT:0>],<EOL>[<NUM_LIT:100> - b for b in list(zip(*b_values))[<NUM_LIT:2>]], c='<STR_LIT:r>')<EOL>ax1.set_ylabel('<STR_LIT>')<EOL>plt.xlabel('<STR_LIT>')<EOL>ax2 = ax1.twinx()<EOL>ax2.set_ylabel('<STR_LIT>')<EOL>n_ev = ax2.scatter(list(zip(*b_values))[<NUM_LIT:0>], list(zip(*b_values))[<NUM_LIT:3>],<EOL>c='<STR_LIT:g>')<EOL>fig.legend((b_vals, resid, n_ev),<EOL>('<STR_LIT>', '<STR_LIT>', '<STR_LIT>'),<EOL>'<STR_LIT>')<EOL>ax1.set_title('<STR_LIT>')<EOL>plt.show()<EOL><DEDENT>return b_values<EOL>
Calculate the b-value for a range of completeness magnitudes.

Calculates a power-law fit to given magnitudes for each completeness
magnitude.  Plots the b-values and residuals for the fitted catalogue
against the completeness values. Computes fits using numpy.polyfit,
which uses a least-squares technique.

:type magnitudes: list
:param magnitudes: Magnitudes to compute the b-value for.
:type completeness: list
:param completeness: list of completeness values to compute b-values for.
:type max_mag: float
:param max_mag: Maximum magnitude to attempt to fit in magnitudes.
:type plotvar: bool
:param plotvar: Turn plotting on or off.

:rtype: list
:return: List of tuples of (completeness, b-value, residual,\
    number of magnitudes used)

.. rubric:: Example

>>> from obspy.clients.fdsn import Client
>>> from obspy import UTCDateTime
>>> from eqcorrscan.utils.mag_calc import calc_b_value
>>> 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.magnitudes[0].mag for event in catalog]
>>> b_values = calc_b_value(magnitudes, completeness=np.arange(3, 7, 0.2),
...                         plotvar=False)
>>> round(b_values[4][1])
1.0
>>> # We can set a maximum magnitude:
>>> b_values = calc_b_value(magnitudes, completeness=np.arange(3, 7, 0.2),
...                         plotvar=False, max_mag=5)
>>> round(b_values[4][1])
1.0