def profile_mask(features, q=2.0, zap_channels=[]):
    """ Compute profile mask, flagging outliers.

    Parameters
    ----------
    features: pandas.DataFrame
        Output of the featurize() function.
    q: float, optional (default: 2.0)
        Parameter that controls the min and max values that define the 'inlier' or
        'normality' range. For every feature, the first and third quartiles (Q1 and Q3)
        are calculated, and R = Q3 - Q1 is the interquartile range. The min and max
        acceptable values are then defined as:

        vmin = Q1 - q x R
        vmax = Q3 + q x R

        The original recommendation of Tukey is q = 1.5.
    zap_channels: list or array of ints
        Frequency channel indices to forcibly flag as outliers. These are excluded
        from the analysis as well, and they will not be considered when calculating
        the inlier range for every feature.

    Returns
    -------
    stats: pandas.DataFrame
        DataFrame containing the results. Columns are individual features, and rows are 
        individual statistics (quartiles, iqr, min and max values).
    mask: ndarray
        A two dimensional boolean numpy array, with shape (num_subints, num_channels).
        A value of True means that the associated profile is an outlier.
    """
    do_zap = zap_channels is not None and len(zap_channels) > 0
    if do_zap:
        X = features.drop(zap_channels, level=1)
    else:
        X = features
    stats = X.quantile([0.25, 0.5, 0.75])
    stats = stats.rename({(0.25): 'q1', (0.5): 'med', (0.75): 'q3'})
    stats.loc['iqr'] = stats.loc['q3'] - stats.loc['q1']
    stats.loc['vmin'] = stats.loc['q1'] - q * stats.loc['iqr']
    stats.loc['vmax'] = stats.loc['q3'] + q * stats.loc['iqr']
    mask = (features > stats.loc['vmax']) | (features < stats.loc['vmin'])
    mask = mask.sum(axis=1).values.astype(bool).reshape(features.index.levshape
        )
    if do_zap:
        mask[:, (zap_channels)] = True
    return stats, mask
