docs: def _calc_inst_pmf(self):
          """Calculate the epsilon-greedy instrumental distribution"""
          # Easy vars
          t = self.t_
          epsilon = self.epsilon
          alpha = self.alpha
          preds = self._preds_avg_in_strata
          weights = self.strata.weights_[:,np.newaxis]
          p1 = self._BB_model.theta_[:,np.newaxis]
          p0 = 1 - p1
          if t==0:
              F = self._F_guess[self.opt_class]
          else:
              F = self._estimate[t - 1, self.opt_class]
              # Fill in non-finite estimates with the initial guess
              nonfinite = ~np.isfinite(F)
              F[nonfinite] = self._F_guess[self.opt_class][nonfinite]
  
          # Calculate optimal instrumental pmf
          sqrt_arg = np.sum(preds * (alpha**2 * F**2 * p0 + (1 - F)**2 * p1) + \
                            (1 - preds) * (1 - alpha)**2 * F**2 * p1, \
                            axis=1, keepdims=True) #: sum is over classifiers
          inst_pmf = weights * np.sqrt(sqrt_arg)
          # Normalize
          inst_pmf /= np.sum(inst_pmf)
          # Epsilon-greedy: (1 - epsilon) q + epsilon * p
          inst_pmf *= (1 - epsilon)
          inst_pmf += epsilon * weights
  
          if self.record_inst_hist:
              self._inst_pmf[:,t] = inst_pmf.ravel()
          else:
              self._inst_pmf = inst_pmf.ravel()