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()