docs: 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)