'Return probability estimates for the test vectors X.
Parameters
X : {array-like, sparse matrix}, shape = [n_samples, n_features]
Input vectors, where n_samples is the number of samples
and n_features is the number of features.
Returns
P : array-like or list of array-lke of shape = [n_samples, n_classes]
Returns the probability of the sample for each class in
the model, where classes are ordered arithmetically, for each
output.'

def predict_proba(self, X):

    check_is_fitted(self, 'classes_')
    X = check_array(X, accept_sparse=['csr', 'csc', 'coo'], force_all_finite=False)
    n_samples = int(X.shape[0])
    rs = check_random_state(self.random_state)
    n_classes_ = self.n_classes_
    classes_ = self.classes_
    class_prior_ = self.class_prior_
    constant = self.constant
    if ((self.n_outputs_ == 1) and (not self.output_2d_)):
        n_classes_ = [n_classes_]
        classes_ = [classes_]
        class_prior_ = [class_prior_]
        constant = [constant]
    P = []
    for k in range(self.n_outputs_):
        if (self.strategy == 'most_frequent'):
            ind = class_prior_[k].argmax()
            out = np.zeros((n_samples, n_classes_[k]), dtype=np.float64)
            out[:, ind] = 1.0
        elif (self.strategy == 'prior'):
            out = (np.ones((n_samples, 1)) * class_prior_[k])
        elif (self.strategy == 'stratified'):
            out = rs.multinomial(1, class_prior_[k], size=n_samples)
        elif (self.strategy == 'uniform'):
            out = np.ones((n_samples, n_classes_[k]), dtype=np.float64)
            out /= n_classes_[k]
        elif (self.strategy == 'constant'):
            ind = np.where((classes_[k] == constant[k]))
            out = np.zeros((n_samples, n_classes_[k]), dtype=np.float64)
            out[:, ind] = 1.0
        P.append(out)
    if ((self.n_outputs_ == 1) and (not self.output_2d_)):
        P = P[0]
    return P
