def overfit(self, tau=None, plot=True, clobber=False, w=9, **kwargs):
        r"""
        Compute the masked & unmasked overfitting metrics for the light curve.

        This routine injects a transit model given by `tau` at every cadence
        in the light curve and recovers the transit depth when (1) leaving
        the transit unmasked and (2) masking the transit prior to performing
        regression.

        :param tau: A function or callable that accepts two arguments, \
               `time` and `t0`, and returns an array corresponding to a \
               zero-mean, unit depth transit model centered at \
               `t0` and evaluated at `time`. \
               The easiest way to provide this is to use an instance of \
               :py:class:`everest.transit.TransitShape`. Default is \
               :py:class:`everest.transit.TransitShape(dur=0.1)`, a transit \
               with solar-like limb darkening and a duratio of 0.1 days.
        :param bool plot: Plot the results as a PDF? Default :py:obj:`True`
        :param bool clobber: Overwrite the results if present? Default \
               :py:obj:`False`
        :param int w: The size of the masking window in cadences for \
               computing the masked overfitting metric. Default `9` \
               (about 4.5 hours for `K2` long cadence).

        :returns: An instance of `everest.basecamp.Overfitting`.
        """
        fname = os.path.join(self.dir, self.name + '_overfit.npz')
        figname = os.path.join(self.dir, self.name)

        # Compute
        if not os.path.exists(fname) or clobber:

            # Baseline
            med = np.nanmedian(self.fraw)

            # Default transit model
            if tau is None:
                tau = TransitShape(dur=0.1)

            # The overfitting metrics
            O1 = [None for brkpt in self.breakpoints]
            O2 = [None for brkpt in self.breakpoints]
            O3 = [None for brkpt in self.breakpoints]
            O4 = [None for brkpt in self.breakpoints]
            O5 = [None for brkpt in self.breakpoints]

            # Loop over all chunks
            for b, brkpt in enumerate(self.breakpoints):

                # Masks for current chunk
                m = self.get_masked_chunk(b, pad=False)
                time = self.time[m]
                ferr = self.fraw_err[m] / med
                y = self.fraw[m] / med - 1

                # The metrics we're computing here
                O1[b] = np.zeros(len(y)) * np.nan
                O2[b] = np.zeros(len(y)) * np.nan
                O3[b] = np.zeros(len(y)) * np.nan
                O4[b] = np.zeros(len(y)) * np.nan
                O5[b] = np.zeros(len(y)) * np.nan

                # Compute the astrophysical covariance and its inverse
                log.info("Computing the covariance...")
                if self.kernel == 'Basic':
                    wh, am, ta = self.kernel_params
                    wh /= med
                    am /= med
                    kernel_params = [wh, am, ta]
                elif self.kernel == 'QuasiPeriodic':
                    wh, am, ga, pe = self.kernel_params
                    wh /= med
                    am /= med
                    kernel_params = [wh, am, ga, pe]
                K = GetCovariance(self.kernel, kernel_params, time, ferr)
                Kinv = cho_solve((cholesky(K), False), np.eye(len(time)))

                # Loop over all orders
                log.info("Computing some large matrices...")
                X = [None for n in range(self.pld_order)]
                XL = [None for n in range(self.pld_order)]
                XLX = [None for n in range(self.pld_order)]
                for n in range(self.pld_order):
                    if (self.lam_idx >= n) and (self.lam[b][n] is not None):
                        X[n] = self.X(n, m, **kwargs)
                        XL[n] = (self.lam[b][n] / med ** 2) * X[n]
                        XLX[n] = np.dot(XL[n], X[n].T)
                X = np.hstack(X)
                XL = np.hstack(XL)
                XLX = np.sum(XLX, axis=0)

                # The full covariance
                C = XLX + K

                # The unmasked linear problem
                log.info("Solving the unmasked linear problem...")
                m = np.dot(XLX, np.linalg.solve(C, y))
                m -= np.nanmedian(m)
                f = y - m
                R = np.linalg.solve(C, XLX.T).T

                # The masked linear problem
                log.info("Solving the masked linear problem...")
                A = MaskSolve(C, y, w=w)

                # Now loop through and compute the metric
                log.info("Computing the overfitting metrics...")
                for n in prange(len(y)):

                    #
                    # *** Unmasked overfitting metric ***
                    #

                    # Evaluate the sparse transit model
                    TAU = tau(time, t0=time[n])
                    i = np.where(TAU < 0)[0]
                    TAU = TAU.reshape(-1, 1)

                    # Fast sparse algebra
                    AA = np.dot(np.dot(TAU[i].T, Kinv[i, :][:, i]), TAU[i])
                    BB = np.dot(TAU[i].T, Kinv[i, :])
                    CC = TAU - np.dot(R[:, i], TAU[i])
                    O1[b][n] = AA
                    O2[b][n] = np.dot(BB, CC)
                    O3[b][n] = np.dot(BB, f)
                    O4[b][n] = np.dot(BB, y)

                    #
                    # *** Masked overfitting metric ***
                    #

                    # The current mask and mask centerpoint
                    mask = np.arange(n, n + w)
                    j = n + (w + 1) // 2 - 1
                    if j >= len(y) - w:
                        continue

                    # The regularized design matrix
                    # This is the same as
                    # XLmX[:, n - 1] = \
                    #   np.dot(XL, np.delete(X, mask, axis=0).T)[:, n - 1]
                    if n == 0:
                        XLmX = np.dot(XL, np.delete(X, mask, axis=0).T)
                    else:
                        XLmX[:, n - 1] = np.dot(XL, X[n - 1, :].T)

                    # The linear solution to this step
                    m = np.dot(XLmX, A[n])

                    # Evaluate the sparse transit model
                    TAU = tau(time, t0=time[j])
                    i = np.where(TAU < 0)[0]
                    TAU = TAU[i].reshape(-1, 1)

                    # Dot the transit model in
                    den = np.dot(np.dot(TAU.T, Kinv[i, :][:, i]), TAU)
                    num = np.dot(TAU.T, Kinv[i, :])

                    # Compute the overfitting metric
                    # Divide this number by a depth
                    # to get the overfitting for that
                    # particular depth.
                    O5[b][j] = -np.dot(num, y - m) / den

            # Save!
            np.savez(fname, O1=O1, O2=O2, O3=O3, O4=O4, O5=O5)

        else:

            data = np.load(fname)
            O1 = data['O1']
            O2 = data['O2']
            O3 = data['O3']
            O4 = data['O4']
            O5 = data['O5']

        # Plot
        if plot and (clobber or not os.path.exists(figname + '_overfit.pdf')):

            log.info("Plotting the overfitting metrics...")

            # Masked time array
            time = self.apply_mask(self.time)

            # Plot the final corrected light curve
            ovr = OVERFIT()
            self.plot_info(ovr)

            # Loop over the two metrics
            for kind, axes, axesh in zip(['unmasked', 'masked'],
                                         [ovr.axes1, ovr.axes2],
                                         [ovr.axes1h, ovr.axes2h]):

                # Loop over three depths
                for depth, ax, axh in zip([0.01, 0.001, 0.0001], axes, axesh):

                    # Get the metric
                    if kind == 'unmasked':
                        metric = 1 - (np.hstack(O2) +
                                      np.hstack(O3) / depth) / np.hstack(O1)
                        color = 'r'
                    elif kind == 'masked':
                        metric = np.hstack(O5) / depth
                        color = 'b'
                    else:
                        raise ValueError("Invalid metric.")

                    # Median and median absolute deviation
                    med = np.nanmedian(metric)
                    mad = np.nanmedian(np.abs(metric - med))

                    # Plot the metric as a function of time
                    ax.plot(time, metric, 'k.', alpha=0.5, ms=2)
                    ax.plot(time, metric, 'k-', alpha=0.1, lw=0.5)
                    ylim = (-0.2, 1.0)
                    ax.margins(0, None)
                    ax.axhline(0, color='k', lw=1, alpha=0.5)
                    ax.set_ylim(*ylim)
                    if kind == 'masked' and depth == 0.0001:
                        ax.set_xlabel('Time (days)', fontsize=14)
                    else:
                        ax.set_xticklabels([])

                    # Plot the histogram
                    rng = (max(ylim[0], np.nanmin(metric)),
                           min(ylim[1], np.nanmax(metric)))
                    axh.hist(metric, bins=30, range=rng,
                             orientation="horizontal",
                             histtype="step", fill=False, color='k')
                    axh.axhline(med, color=color, ls='-', lw=1)
                    axh.axhspan(med - mad, med + mad, color=color, alpha=0.1)
                    axh.axhline(0, color='k', lw=1, alpha=0.5)
                    axh.yaxis.tick_right()
                    axh.set_ylim(*ax.get_ylim())
                    axh.set_xticklabels([])

                    bbox = dict(fc="w", ec="1", alpha=0.5)
                    info = r"$\mathrm{med}=%.3f$" % med + \
                        "\n" + r"$\mathrm{mad}=%.3f$" % mad
                    axh.annotate(info, xy=(0.1, 0.925),
                                 xycoords='axes fraction',
                                 ha="left", va="top", bbox=bbox, color=color)

                    bbox = dict(fc="w", ec="1", alpha=0.95)
                    ax.annotate("%s overfitting metric" % kind,
                                xy=(1-0.035, 0.92),
                                xycoords='axes fraction',
                                ha='right', va='top',
                                bbox=bbox, color=color)

            pl.figtext(0.025, 0.77, "depth = 0.01", rotation=90,
                       ha='left', va='center', fontsize=18)
            pl.figtext(0.025, 0.48, "depth = 0.001", rotation=90,
                       ha='left', va='center', fontsize=18)
            pl.figtext(0.025, 0.19, "depth = 0.0001", rotation=90,
                       ha='left', va='center', fontsize=18)

            ovr.fig.savefig(figname + '_overfit.pdf')
            log.info("Saved plot to %s_overfit.pdf" % figname)
            pl.close()

        return Overfitting(O1, O2, O3, O4, O5, figname + '_overfit.pdf')