docs: def fit_fb_calibration(df, calibration):
      '''
      Fit feedback calibration data to solve for values of `C_fb[:]` and
      `R_fb[:]`.
  
      Returns a `pandas.DataFrame` indexed by the feedback resistor/capacitance
      index, and with the following columns:
       - Model: Either with parasitic capacitance term or not.
       - N: Number of samples used for fit.
       - F: F-value
       - p-value: p-value from Chi squared test.
       - R_fb: Feedback resistor value based on fit.
       - R-CI %: Confidence interval for feedback resistor value.
       - C_fb: Feedback capacitor value based on fit (0 if no-capacitance model
         is used).
       - C-CI %: Confidence interval for feedback capacitance value.
  
      __N.B.__ This function does not actually _update_ the calibration, it only
      performs the fit.
      See `apply_calibration`.
      '''
      # Set initial guesses for the feedback parameters.
      R_fb = pd.Series([2e2, 2e3, 2e4, 2e5, 2e6])
      C_fb = pd.Series(len(calibration.C_fb) * [50e-12])
  
      # Error function.
      def error(p0, df, calibration):
          # Impedance of the reference resistor on the HV attenuator circuit.
          Z = 10e6
          R_fb = p0[0]
  
          # If the parameter vector only contains one variable, the capacitance
          # is zero
          if len(p0) == 2:
              C_fb = p0[1]
          else:
              C_fb = 0
  
          R_hv = calibration.R_hv[df.hv_resistor.values]
          C_hv = calibration.C_hv[df.hv_resistor.values]
  
          # Solve feedback transfer function for the actuation voltage, _(i.e.,
          # `V1`)_, based on the high-voltage measurements.
          # Note that the transfer function definition depends on the hardware
          # version.
          V_actuation = compute_from_transfer_function(calibration.hw_version
                                                       .major, 'V1', V2=df.V_hv,
                                                       R1=Z, R2=R_hv, C2=C_hv,
                                                       f=df.frequency)
  
          # Solve feedback transfer function for the expected impedance feedback
          # voltage, _(i.e., `V2`)_, based on the actuation voltage, the proposed
          # values for `R2` and `C2`, and the reported `C1` value from the
          # feedback measurements.
          # Note that the transfer function definition depends on the hardware
          # version.
          # __NB__ If we do not specify a value for `R1`, a symbolic value of
          # infinity is used.  However, in this case, we have `R1` in both the
          # numerator and denominator.  The result is a value of zero returned
          # regardless of the values of the other arguments.  We avoid this issue
          # by specifying a *very large* value for `R1`.
          # TODO Update comment if this works...
          V_impedance = compute_from_transfer_function(calibration.hw_version
                                                       .major, 'V2',
                                                       V1=V_actuation,
                                                       C1=df.test_capacitor,
                                                       R2=R_fb, C2=C_fb,
                                                       f=df.frequency)
          return df.V_fb - V_impedance
  
      # Perform a nonlinear least-squares fit of the data.
      def fit_model(p0, df, calibration):
          p1, cov_x, infodict, mesg, ier = scipy.optimize.leastsq(
              error, p0, args=(df, calibration), full_output=True)
          p1 = np.abs(p1)
          E = error(p1, df, calibration)
          return p1, E, cov_x
  
      CI = []
  
      feedback_records = []
      # Fit feedback parameters for each feedback resistor.
      for i in range(len(calibration.R_fb)):
          # Only include data points for the given feedback resistor (and where
          # `hv_resistor` is a valid index).
          df_i = df.loc[(df.fb_resistor == i)].dropna()
  
          if df_i.shape[0] < 2:
              CI.append([0, 0])
              continue
  
          # Fit the data assuming no parasitic capacitance (model 1).
          p0_1 = [R_fb[i]]
          p1_1, E_1, cov_x_1 = fit_model(p0_1, df_i, calibration)
          df_1 = (len(E_1) - len(p0_1))
          chi2_1 = np.sum(E_1 ** 2)
          chi2r_1 = chi2_1 / (df_1 - 1)
  
          # fit the data including parasitic capacitance (model 2)
          p0_2 = [R_fb[i], C_fb[i]]
          p1_2, E_2, cov_x_2 = fit_model(p0_2, df_i, calibration)
          df_2 = (len(E_2) - len(p0_2))
          chi2_2 = np.sum(E_2 ** 2)
          chi2r_2 = chi2_2 / (df_2 - 1)
  
          # do an F-test to compare the models
          F = (chi2_1 - chi2_2) / chi2r_2
          p_value = scipy.stats.f.cdf(F, 1, df_2-1)
  
          # if the p_value is > 0.95, we assume that the capacitive term is
          # necessary
          if p_value > .95 and cov_x_2 is not None:
              model = 'w/Parasitic C'
              chi2r = chi2r_2
              R_fb_i = p1_2[0]
              C_fb_i = p1_2[1]
              CI.append((100 * np.sqrt(chi2r_2 * np.diag(cov_x_2)) / p1_2))
          else:  # otherwise, set the capacitance to zero
              model = 'w/o Parasitic C'
              chi2r = chi2r_2
              R_fb_i = p1_1[0]
              C_fb_i = 0
              if cov_x_1 is None:
                  cov_x_1 = [0]
              CI.append((100 * np.sqrt(chi2r_1 * np.diag(cov_x_1)) /
                         p1_1).tolist() + [0])
          feedback_records.append([int(i), model, df_i.shape[0], R_fb_i, CI[i][0],
                                   C_fb_i, CI[i][1], F, (1e3 * np.sqrt(chi2r)),
                                   p_value])
  
      calibration_df = pd.DataFrame(feedback_records,
                                    columns=['fb_resistor', 'Model', 'N', 'R_fb', 'R-CI %',
                                             'C_fb', 'C-CI %', 'F',
                                             'sqrt(Chi2r*sigma^2)', 'p-value'])
      return calibration_df