def Statistics(season=None, clobber=False, model='nPLD', injection=False,
               compare_to='kepler', plot=True, cadence='lc', planets=False,
               **kwargs):
    '''
    Computes and plots the CDPP statistics comparison between `model`
    and `compare_to` for all long cadence light curves in a given campaign

    :param season: The campaign number or list of campaign numbers. \
           Default is to plot all campaigns
    :param bool clobber: Overwrite existing files? Default :py:obj:`False`
    :param str model: The :py:obj:`everest` model name
    :param str compare_to: The :py:obj:`everest` model name or other \
           K2 pipeline name
    :param bool plot: Default :py:obj:`True`
    :param bool injection: Statistics for injection tests? Default \
           :py:obj:`False`
    :param bool planets: Statistics for known K2 planets? \
           Default :py:obj:`False`

    '''

    # Multi-mission compatibility
    campaign = season

    # Is this short cadence?
    if cadence == 'sc':
        return ShortCadenceStatistics(campaign=campaign, clobber=clobber,
                                      model=model, plot=plot, **kwargs)

    # Check the campaign
    if campaign is None:
        campaign = 0

    # Planet hosts only?
    if planets:
        return PlanetStatistics(model=model, compare_to=compare_to, **kwargs)

    # Is this an injection run?
    if injection:
        return InjectionStatistics(campaign=campaign, clobber=clobber,
                                   model=model, plot=plot, **kwargs)

    # Compute the statistics
    sub = np.array([s[0] for s in GetK2Campaign(campaign)], dtype=int)
    outfile = os.path.join(EVEREST_SRC, 'missions', 'k2',
                           'tables', 'c%02d_%s.cdpp' % (int(campaign), model))
    if clobber or not os.path.exists(outfile):
        with open(outfile, 'w') as f:
            print("EPIC               Kp           Raw CDPP     Everest CDPP" +
                  "      Validation        Outliers[1]     Outliers[2]     " +
                  "Datapoints     Saturated", file=f)
            print("---------          ------       ---------    ------------" +
                  "      ----------        -----------     -----------     " +
                  "----------     ---------", file=f)
            all = GetK2Campaign(int(campaign))
            stars = np.array([s[0] for s in all], dtype=int)
            kpmgs = np.array([s[1] for s in all], dtype=float)
            for i, _ in enumerate(stars):
                sys.stdout.write('\rProcessing target %d/%d...' %
                                 (i + 1, len(stars)))
                sys.stdout.flush()
                nf = os.path.join(EVEREST_DAT, 'k2', 'c%02d' % campaign,
                                  ('%09d' % stars[i])[:4] + '00000',
                                  ('%09d' % stars[i])[4:], model + '.npz')
                try:
                    data = np.load(nf)

                    # Remove NaNs and flagged cadences
                    flux = np.delete(data['fraw'] - data['model'], np.array(
                        list(set(np.concatenate([data['nanmask'],
                                                 data['badmask']])))))
                    # Iterative sigma clipping to get 5 sigma outliers
                    inds = np.array([], dtype=int)
                    m = 1
                    while len(inds) < m:
                        m = len(inds)
                        ff = SavGol(np.delete(flux, inds))
                        med = np.nanmedian(ff)
                        MAD = 1.4826 * np.nanmedian(np.abs(ff - med))
                        inds = np.append(inds, np.where(
                            (ff > med + 5. * MAD) | (ff < med - 5. * MAD))[0])
                    nout = len(inds)
                    ntot = len(flux)

                    # HACK: Backwards compatibility fix
                    try:
                        cdpp = data['cdpp'][()]
                    except KeyError:
                        cdpp = data['cdpp6'][()]

                    print("{:>09d} {:>15.3f} {:>15.3f} {:>15.3f} {:>15.3f} {:>15d} {:>15d} {:>15d} {:>15d}".format(
                          stars[i], kpmgs[i], data['cdppr'][()], cdpp,
                          data['cdppv'][()], len(data['outmask']), nout,
                          ntot, int(data['saturated'])), file=f)
                except:
                    print("{:>09d} {:>15.3f} {:>15.3f} {:>15.3f} {:>15.3f} {:>15d} {:>15d} {:>15d} {:>15d}".format(
                          stars[i], kpmgs[i], np.nan, np.nan,
                          np.nan, 0, 0, 0, 0), file=f)
            print("")

    if plot:

        # Load all stars
        epic, kp, cdpp6r, cdpp6, cdpp6v, _, out, tot, saturated = np.loadtxt(
            outfile, unpack=True, skiprows=2)
        epic = np.array(epic, dtype=int)
        out = np.array(out, dtype=int)
        tot = np.array(tot, dtype=int)
        saturated = np.array(saturated, dtype=int)

        # Get only stars in this subcampaign
        inds = np.array([e in sub for e in epic])
        epic = epic[inds]
        kp = kp[inds]
        # HACK: Campaign 0 magnitudes are reported only to the nearest tenth,
        # so let's add a little noise to spread them out for nicer plotting
        kp = kp + 0.1 * (0.5 - np.random.random(len(kp)))
        cdpp6r = cdpp6r[inds]
        cdpp6 = cdpp6[inds]
        cdpp6v = cdpp6v[inds]
        out = out[inds]
        tot = tot[inds]
        saturated = saturated[inds]
        sat = np.where(saturated == 1)
        unsat = np.where(saturated == 0)
        if not np.any([not np.isnan(x) for x in cdpp6]):
            raise Exception("No targets to plot.")

        # Control transparency
        alpha_kepler = 0.03
        alpha_unsat = min(0.1, 2000. / (1 + len(unsat[0])))
        alpha_sat = min(1., 180. / (1 + len(sat[0])))

        # Get the comparison model stats
        if compare_to.lower() == 'everest1':
            epic_1, cdpp6_1 = np.loadtxt(
                              os.path.join(EVEREST_SRC, 'missions', 'k2',
                                           'tables', 'c%02d_everest1.cdpp'
                                           % int(campaign)), unpack=True)
            cdpp6_1 = sort_like(cdpp6_1, epic, epic_1)
            # Outliers
            epic_1, out_1, tot_1 = np.loadtxt(
                                   os.path.join(EVEREST_SRC, 'missions', 'k2',
                                                'tables', 'c%02d_everest1.out'
                                                % int(campaign)), unpack=True)
            out_1 = sort_like(out_1, epic, epic_1)
            tot_1 = sort_like(tot_1, epic, epic_1)
        elif compare_to.lower() == 'k2sc':
            epic_1, cdpp6_1 = np.loadtxt(
                              os.path.join(EVEREST_SRC, 'missions', 'k2',
                                           'tables', 'c%02d_k2sc.cdpp' %
                                           int(campaign)), unpack=True)
            cdpp6_1 = sort_like(cdpp6_1, epic, epic_1)
            # Outliers
            epic_1, out_1, tot_1 = np.loadtxt(
                                   os.path.join(EVEREST_SRC, 'missions', 'k2',
                                                'tables', 'c%02d_k2sc.out'
                                                % int(campaign)), unpack=True)
            out_1 = sort_like(out_1, epic, epic_1)
            tot_1 = sort_like(tot_1, epic, epic_1)
        elif compare_to.lower() == 'k2sff':
            epic_1, cdpp6_1 = np.loadtxt(
                              os.path.join(EVEREST_SRC, 'missions', 'k2',
                                           'tables', 'c%02d_k2sff.cdpp'
                                           % int(campaign)), unpack=True)
            cdpp6_1 = sort_like(cdpp6_1, epic, epic_1)
            # Outliers
            epic_1, out_1, tot_1 = np.loadtxt(
                                   os.path.join(EVEREST_SRC, 'missions', 'k2',
                                                'tables', 'c%02d_k2sff.out'
                                                % int(campaign)), unpack=True)
            out_1 = sort_like(out_1, epic, epic_1)
            tot_1 = sort_like(tot_1, epic, epic_1)
        elif compare_to.lower() == 'kepler':
            kic, kepler_kp, kepler_cdpp6 = np.loadtxt(
                                           os.path.join(EVEREST_SRC,
                                                        'missions', 'k2',
                                                        'tables',
                                                        'kepler.cdpp'),
                                           unpack=True)
        else:
            compfile = os.path.join(EVEREST_SRC, 'missions', 'k2',
                                    'tables', 'c%02d_%s.cdpp' %
                                    (int(campaign), compare_to))

            epic_1, _, _, cdpp6_1, _, _, out_1, tot_1, saturated = np.loadtxt(
                compfile, unpack=True, skiprows=2)
            epic_1 = np.array(epic_1, dtype=int)
            inds = np.array([e in sub for e in epic_1])
            epic_1 = epic_1[inds]
            cdpp6_1 = cdpp6_1[inds]
            out_1 = out_1[inds]
            tot_1 = tot_1[inds]
            cdpp6_1 = sort_like(cdpp6_1, epic, epic_1)
            out_1 = sort_like(out_1, epic, epic_1)
            tot_1 = sort_like(tot_1, epic, epic_1)

        # ------ 1. Plot cdpp vs. mag
        if compare_to.lower() != 'kepler':
            fig = pl.figure(figsize=(16, 5))
            ax = [pl.subplot2grid((120, 120), (0,  0), colspan=35,
                                  rowspan=120),
                  pl.subplot2grid((120, 120), (0,  40),
                                  colspan=35, rowspan=120),
                  pl.subplot2grid((120, 120), (0,  80),
                                  colspan=35, rowspan=55),
                  pl.subplot2grid((120, 120), (65,  80), colspan=35,
                                  rowspan=55)]
        else:
            fig = pl.figure(figsize=(12, 5))
            ax = [pl.subplot2grid((120, 75), (0,  0), colspan=35, rowspan=120),
                  None,
                  pl.subplot2grid((120, 75), (0,  40), colspan=35, rowspan=55),
                  pl.subplot2grid((120, 75), (65,  40), colspan=35,
                                  rowspan=55)]
        fig.canvas.set_window_title(
            'K2 Campaign %s: %s versus %s' % (campaign, model, compare_to))
        fig.subplots_adjust(left=0.05, right=0.95, bottom=0.125, top=0.9)
        bins = np.arange(7.5, 18.5, 0.5)
        if compare_to.lower() != 'kepler':
            ax[0].scatter(kp[unsat], cdpp6_1[unsat], color='y',
                          marker='.', alpha=alpha_unsat)
            ax[0].scatter(kp[sat], cdpp6_1[sat], color='y',
                          marker='s', alpha=alpha_sat, s=5)
            ax[0].scatter(kp[unsat], cdpp6[unsat], color='b',
                          marker='.', alpha=alpha_unsat, picker=True)
            ax[0].scatter(kp[sat], cdpp6[sat], color='b',
                          marker='s', alpha=alpha_sat, s=5, picker=True)
            for y, style in zip([cdpp6_1, cdpp6], ['yo', 'bo']):
                by = np.zeros_like(bins) * np.nan
                for b, bin in enumerate(bins):
                    i = np.where((y > -np.inf) & (y < np.inf) &
                                 (kp >= bin - 0.5) & (kp < bin + 0.5))[0]
                    if len(i) > 10:
                        by[b] = np.median(y[i])
                ax[0].plot(bins, by, style, markeredgecolor='w')
        else:
            ax[0].scatter(kepler_kp, kepler_cdpp6, color='y',
                          marker='.', alpha=alpha_kepler)
            ax[0].scatter(kp, cdpp6, color='b', marker='.',
                          alpha=alpha_unsat, picker=True)
            for x, y, style in zip([kepler_kp, kp], [kepler_cdpp6, cdpp6],
                                   ['yo', 'bo']):
                by = np.zeros_like(bins) * np.nan
                for b, bin in enumerate(bins):
                    i = np.where((y > -np.inf) & (y < np.inf) &
                                 (x >= bin - 0.5) & (x < bin + 0.5))[0]
                    if len(i) > 10:
                        by[b] = np.median(y[i])
                ax[0].plot(bins, by, style, markeredgecolor='w')
        ax[0].set_ylim(-10, 500)
        ax[0].set_xlim(8, 18)
        ax[0].set_xlabel('Kepler Magnitude', fontsize=18)
        ax[0].set_title('CDPP6 (ppm)', fontsize=18)

        # ------ 2. Plot the equivalent of the Aigrain+16 figure
        if compare_to.lower() != 'kepler':
            x = kp
            y = (cdpp6 - cdpp6_1) / cdpp6_1
            yv = (cdpp6v - cdpp6_1) / cdpp6_1
            ax[1].scatter(x[unsat], y[unsat], color='b', marker='.',
                          alpha=alpha_unsat, zorder=-1, picker=True)
            ax[1].scatter(x[sat], y[sat], color='r', marker='.',
                          alpha=alpha_sat, zorder=-1, picker=True)
            ax[1].set_ylim(-1, 1)
            ax[1].set_xlim(8, 18)
            ax[1].axhline(0, color='gray', lw=2, zorder=-99, alpha=0.5)
            ax[1].axhline(0.5, color='gray', ls='--',
                          lw=2, zorder=-99, alpha=0.5)
            ax[1].axhline(-0.5, color='gray', ls='--',
                          lw=2, zorder=-99, alpha=0.5)
            bins = np.arange(7.5, 18.5, 0.5)
            # Bin the CDPP
            by = np.zeros_like(bins) * np.nan
            for b, bin in enumerate(bins):
                i = np.where((y > -np.inf) & (y < np.inf) &
                             (x >= bin - 0.5) & (x < bin + 0.5))[0]
                if len(i) > 10:
                    by[b] = np.median(y[i])
            ax[1].plot(bins[:9], by[:9], 'k--', lw=2)
            ax[1].plot(bins[8:], by[8:], 'k-', lw=2)
            ax[1].set_title(r'Relative CDPP', fontsize=18)
            ax[1].set_xlabel('Kepler Magnitude', fontsize=18)

        # ------ 3. Plot the outliers
        i = np.argsort(out)
        a = int(0.95 * len(out))
        omax = out[i][a]
        if compare_to.lower() != 'kepler':
            j = np.argsort(out_1)
            b = int(0.95 * len(out_1))
            omax = max(omax, out_1[j][b])
        ax[2].hist(out, 25, range=(0, omax), histtype='step', color='b')
        if compare_to.lower() != 'kepler':
            ax[2].hist(out_1, 25, range=(0, omax), histtype='step', color='y')
        ax[2].margins(0, None)
        ax[2].set_title('Number of Outliers', fontsize=18)

        # Plot the total number of data points
        i = np.argsort(tot)
        a = int(0.05 * len(tot))
        b = int(0.95 * len(tot))
        tmin = tot[i][a]
        tmax = tot[i][b]
        if compare_to.lower() != 'kepler':
            j = np.argsort(tot_1)
            c = int(0.05 * len(tot_1))
            d = int(0.95 * len(tot_1))
            tmin = min(tmin, tot_1[j][c])
            tmax = max(tmax, tot_1[j][d])
        ax[3].hist(tot, 25, range=(tmin, tmax), histtype='step', color='b')
        if compare_to.lower() != 'kepler':
            ax[3].hist(tot_1, 25, range=(tmin, tmax),
                       histtype='step', color='y')
        ax[3].margins(0, None)
        ax[3].set_xlabel('Number of Data Points', fontsize=18)

        # Pickable points
        Picker = StatsPicker([ax[0], ax[1]], [kp, kp], [
                             cdpp6, y], epic, model=model,
                             compare_to=compare_to, campaign=campaign)
        fig.canvas.mpl_connect('pick_event', Picker)

        # Show
        pl.show()