'Return probability estimates for the test vector X.
Parameters
X : array-like, shape = (n_samples, n_features)
Returns
C : array-like, shape = (n_samples, n_classes)
Returns the probability of the samples for each class in
the model. The columns correspond to the classes in sorted
order, as they appear in the attribute ``classes_``.'

def predict_proba(self, X):

    check_is_fitted(self, ['X_train_', 'y_train_', 'pi_', 'W_sr_', 'L_'])
    K_star = self.kernel_(self.X_train_, X)
    f_star = K_star.T.dot((self.y_train_ - self.pi_))
    v = solve(self.L_, (self.W_sr_[:, np.newaxis] * K_star))
    var_f_star = (self.kernel_.diag(X) - np.einsum('ij,ij->j', v, v))
    alpha = (1 / (2 * var_f_star))
    gamma = (LAMBDAS * f_star)
    integrals = ((np.sqrt((np.pi / alpha)) * erf((gamma * np.sqrt((alpha / (alpha + (LAMBDAS ** 2))))))) / (2 * np.sqrt(((var_f_star * 2) * np.pi))))
    pi_star = ((COEFS * integrals).sum(axis=0) + (0.5 * COEFS.sum()))
    return np.vstack(((1 - pi_star), pi_star)).T
