def pairwise_reproducibility(df, plot=False):
    """
    Calculate the reproducibility of LA-ICPMS based on unique pairs of repeat analyses.
    
    Pairwise differences are fit with a half-Cauchy distribution, and the median and 
    95% confidence limits are returned for each analyte.
    
    Parameters
    ----------
    df : pandas.DataFrame
        A dataset
    
    plot : bool
        Whether or not to plot the resulting error distributions.
    
    Returns
    -------
    pdiffs : pandas.DataFrame
        Unique pairwise differences for all analytes.
    rep_dists : dict of scipy.stats.halfcauchy
        Half-Cauchy distribution objects fitted to the
        differences.
    rep_stats : dict of tuples
        The 50% and 95% quantiles of the half-cauchy
        distribution.
    (fig, axs) : matplotlib objects
        The figure. If not made, returnes (None, None) placeholder
    
    """
    
    ans = df.columns.values
    pdifs = []
    
    # calculate differences between unique pairs
    for ind, d in df.groupby(level=0):
        d.index = d.index.droplevel(0)

        difs = []
        for i, r in d.iterrows():
            t = d.loc[i+1:, :]
            difs.append(t[ans] - r[ans])

        pdifs.append(pd.concat(difs))
    pdifs = pd.concat(pdifs).abs()

    # calculate stats
    rep_stats = {}
    rep_dists = {}
    errfn = stats.halfcauchy
    
    for a in ans:
        d = pdifs.loc[:, a].dropna().values
        hdist = errfn.fit(d, floc=0)
        rep_dists[a] = errfn(*hdist)
        rep_stats[a] = rep_dists[a].ppf((0.5, 0.95))
    
    # make plot
    if not plot:
        return pdifs, rep_dists, rep_stats, (None, None)
    
    fig, axs = plt.subplots(1, len(ans), figsize=[len(ans) * 2, 2])
    for a, ax in zip(ans, axs):
        d = pdifs.loc[:, a].dropna().values
        hist, edges, _ = ax.hist(d, 30)
        ax.plot(edges, rep_dists[a].pdf(edges) * (sum(hist) * np.mean(np.diff(edges))))
        ax.set_title(a, loc='left')

    return pdifs, rep_dists, rep_stats, (fig, axs)
Calculate the reproducibility of LA-ICPMS based on unique pairs of repeat analyses.
    
    Pairwise differences are fit with a half-Cauchy distribution, and the median and 
    95% confidence limits are returned for each analyte.
    
    Parameters
    ----------
    df : pandas.DataFrame
        A dataset
    
    plot : bool
        Whether or not to plot the resulting error distributions.
    
    Returns
    -------
    pdiffs : pandas.DataFrame
        Unique pairwise differences for all analytes.
    rep_dists : dict of scipy.stats.halfcauchy
        Half-Cauchy distribution objects fitted to the
        differences.
    rep_stats : dict of tuples
        The 50% and 95% quantiles of the half-cauchy
        distribution.
    (fig, axs) : matplotlib objects
        The figure. If not made, returnes (None, None) placeholder