def _initialise_pmf(self):
        """Calculate the epsilon-greedy instrumental distribution"""
        # Easy vars
        epsilon = self.epsilon
        alpha = self.alpha
        preds = self.predictions
        p1 = self._probs_avg_opt_class
        p0 = 1 - p1
        n_items = self._n_items
        F = self._F_guess

        # 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) #: sum is over classifiers
        self._inst_pmf = np.sqrt(sqrt_arg)
        # Normalize
        self._inst_pmf /= np.sum(self._inst_pmf)
        # Epsilon-greedy: (1 - epsilon) q + epsilon * p
        self._inst_pmf *= (1 - epsilon)
        self._inst_pmf += epsilon * 1/n_items