Plot histograms of data residuals and data error weighting

        TODO:
            * add percentage of data below/above the RMS value

def plot_eps_data_hist(self, dfs):
        """Plot histograms of data residuals and data error weighting

        TODO:
            * add percentage of data below/above the RMS value
        """
        # check if this is a DC inversion
        if 'datum' in dfs[0]:
            dc_inv = True
        else:
            dc_inv = False

        nr_y = len(dfs)
        size_y = 5 / 2.54 * nr_y
        if dc_inv:
            nr_x = 1
        else:
            nr_x = 3
        size_x = 15 / 2.54

        fig, axes = plt.subplots(nr_y, nr_x, figsize=(size_x, size_y))
        axes = np.atleast_2d(axes)

        # plot initial data errors
        df = dfs[0]
        if dc_inv:
            ax = axes[0, 0]
            ax.hist(
                df['datum'] / df['eps_r'],
                100,
            )
            ax.set_xlabel(r'$-log(|R|) / \epsilon_r$')
            ax.set_ylabel(r'count')
        else:
            # complex inversion
            ax = axes[0, 0]
            ax.hist(
                df['-log(|R|)'] / df['eps'],
                100,
            )
            ax.set_xlabel(r'$-log(|R|)$')
            ax.set_ylabel(r'count')

            ax = axes[0, 1]
            ax.hist(
                df['-log(|R|)'] / df['eps_r'],
                100,
            )
            ax.set_xlabel(r'$-log(|R|) / \epsilon_r$')
            ax.set_ylabel(r'count')

            ax = axes[0, 2]
            phase_data = df['-Phase(rad)'] / df['eps_p']
            if not np.all(np.isinf(phase_data) | np.isnan(phase_data)):
                ax.hist(
                    phase_data,
                    100,
                )
            ax.set_xlabel(r'$-\phi[rad] / \epsilon_p$')
            ax.set_ylabel(r'count')

        # iterations
        for it, df in enumerate(dfs[1:]):
            ax = axes[1 + it, 0]
            ax.hist(
                df['psi'],
                100
            )
            rms = np.sqrt(
                1 / df['psi'].shape[0] *
                np.sum(
                    df['psi'] ** 2
                )
            )
            ax.axvline(rms, color='k', linestyle='dashed')
            ax.set_title('iteration: {0}'.format(it))
            ax.set_xlabel('psi')
            ax.set_ylabel(r'count')

            ax = axes[1 + it, 1]
            Rdat = df['Re(d)']
            Rmod = df['Re(f(m))']

            ax.scatter(
                Rdat,
                Rmod,
            )
            ax.set_xlabel(r'$log(R_{data}~[\Omega])$')
            ax.set_ylabel(r'$log(R_{mod}~[\Omega])$')

            ax = axes[1 + it, 2]
            phidat = df['Im(d)']
            phimod = df['Im(f(m))']

            ax.scatter(
                phidat,
                phimod,
            )
            ax.set_xlabel(r'$\phi_{data}~[mrad]$')
            ax.set_ylabel(r'$\phi_{mod}~[mrad]$')

        fig.tight_layout()
        fig.savefig('eps_plot_hist.png', dpi=300)