Lets say my multinomial logistic regression predict that a chance of a sample belonging to a each class is A=0.6, B=0.3, C=0.1 How do I threshold this values to get just binary prediction of a sample belonging to a class, taking in to an account imbalances of classes. I know what I would do if it's just a binary decision (threshold based on classes prevalence), or if the classes are balanced (classify to a class with highest probability). My end goal is to get 3x3 confusion matrix

  • $\begingroup$ I'm a bit unclear on what you are trying to achieve. If you use three different thresholds $T_A, T_B, T_C$ and classify a sample as $A$ if its predicted probability satisfies $p_A\geq T_A$, then you may well end up in situations where both $p_A\geq T_B$ and $p_B\geq T_B$, and what will you do then? I'd recommend classifying everything based on the highest probability, then write up the confusion matrix, then assign costs to each kind of misclassification, look which misclassifications are the main drivers of total costs, finally tune the classifier itself on these misclassifications. $\endgroup$ Sep 29, 2016 at 10:01
  • $\begingroup$ That's exactly my problem. IF I classify everything based on highest probability, I can end up classifying everything as majority class. I don't want to tune the classifier, just the predictions. $\endgroup$
    – rep_ho
    Sep 29, 2016 at 10:04
  • $\begingroup$ Let's go back to binary classification, which outputs a single probability $p$ for membership in the "target" class. If you choose a threshold $T$ and classify a sample as "target" whenever $p\geq T$, then choosing $T$ is part of tuning the classifier - after all, the desired output of a classifier is a classification, not some probability value that you still need to compare to a threshold. Plus: if you want to work with thresholds, how do you plan on dealing with $p_A\geq T_A$ & $p_B\geq T_B$ situations? $\endgroup$ Sep 29, 2016 at 10:09
  • $\begingroup$ let's say in a binary classifier my classes are balanced: A=0.1, B=0.9 and my classifier predicts A=0.4 and B=0.6. Although B has higher probability, I would classify it as A because the predicted probability is higher than a prior. You can tune your threshold to be better with respect to some measure, but I am willing to use this simplified untuned version. I don't know what to do in pA≥TA & pB≥TB situations, that's part of the problem and why I asked the question $\endgroup$
    – rep_ho
    Sep 29, 2016 at 10:21
  • $\begingroup$ Why is your criterion whether or not the predicted probability is higher than a prior? I'd argue that you should set your threshold based on the costs of misclassifications. And yes, classifying everything in the most common class may well be your cost-optimal classifier. (Store employees routinely classify everyone that enters as "honest", unless their "classifier" gets some very specific information.) That's why I say that your threshold (which is part of the classifier) needs to take costs into account. And similarly, costs should influence on how you train your multinomial classifier. $\endgroup$ Sep 29, 2016 at 10:29

2 Answers 2


According to @cangrejo's answer: https://stats.stackexchange.com/a/310956/194535, suppose the original output probability of your model is the vector $v$, and then you can define the prior distribution:

$\pi=(\frac{1}{\theta_1}, \frac{1}{\theta_2},..., \frac{1}{\theta_N})$, for $\theta_i \in (0,1)$ and $\sum_i\theta_i = 1$, where $N$ is the total number of labeled classes, $i$ is the class index.

Take $v' = v \odot \pi$ as the new output probability of your model, where $\odot$ denotes an element-wise product.

Now, your question can be reformulate to this: Finding the $\pi$ which optimize the metrics you have specified (eg. roc_auc_score) from the new output probability model. Once you find it, the $\theta s (\theta_1, \theta_2, ..., \theta_N)$ is your optimal threshold for each classes.

The Code part:

  1. Create a proxyModel class which takes your original model object as an argument and return a proxyModel object. When you called predict_proba() through the proxyModel object, it will calculate new probability automatically based on the threshold you specified:

    class proxyModel():
        def __init__(self, origin_model):
            self.origin_model = origin_model
        def predict_proba(self, x, threshold_list=None):
            # get origin probability
            ori_proba = self.origin_model.predict_proba(x)
            # set default threshold
            if threshold_list is None:
                threshold_list = np.full(ori_proba[0].shape, 1)
            # get the output shape of threshold_list
            output_shape = np.array(threshold_list).shape
            # element-wise divide by the threshold of each classes
            new_proba = np.divide(ori_proba, threshold_list)
            # calculate the norm (sum of new probability of each classes)
            norm = np.linalg.norm(new_proba, ord=1, axis=1)
            # reshape the norm
            norm = np.broadcast_to(np.array([norm]).T, (norm.shape[0],output_shape[0]))
            # renormalize the new probability
            new_proba = np.divide(new_proba, norm)
            return new_proba
        def predict(self, x, threshold_list=None):
            return np.argmax(self.predict_proba(x, threshold_list), axis=1)
  2. Implement a score function:

    def scoreFunc(model, X, y_true, threshold_list):
        y_pred = model.predict(X, threshold_list=threshold_list)
        y_pred_proba = model.predict_proba(X, threshold_list=threshold_list)
        ###### metrics ######
        from sklearn.metrics import accuracy_score
        from sklearn.metrics import roc_auc_score
        from sklearn.metrics import average_precision_score
        from sklearn.metrics import f1_score
        accuracy = accuracy_score(y_true, y_pred)
        roc_auc = roc_auc_score(y_true, y_pred_proba, average='macro')
        pr_auc = average_precision_score(y_true, y_pred_proba, average='macro')
        f1_value = f1_score(y_true, y_pred, average='macro')
        return accuracy, roc_auc, pr_auc, f1_value
  3. Define weighted_score_with_threshold() function, which takes the threshold as input and return weighted score:

    def weighted_score_with_threshold(threshold, model, X_test, Y_test, metrics='accuracy', delta=5e-5):
        # if the sum of thresholds were not between 1+delta and 1-delta, 
        # return infinity (just for reduce the search space of the minimizaiton algorithm, 
        # because the sum of thresholds should be as close to 1 as possible).
        threshold_sum = np.sum(threshold)
        if threshold_sum > 1+delta:
            return np.inf
        if threshold_sum < 1-delta:
            return np.inf
        # to avoid objective function jump into nan solution
        if np.isnan(threshold_sum):
            print("threshold_sum is nan")
            return np.inf
        # renormalize: the sum of threshold should be 1
        normalized_threshold = threshold/threshold_sum
        # calculate scores based on thresholds
        # suppose it'll return 4 scores in a tuple: (accuracy, roc_auc, pr_auc, f1)
        scores = scoreFunc(model, X_test, Y_test, threshold_list=normalized_threshold)    
        scores = np.array(scores)
        weight = np.array([1,1,1,1])
        # Give the metric you want to maximize a bigger weight:
        if metrics == 'accuracy':
            weight = np.array([10,1,1,1])
        elif metrics == 'roc_auc':
            weight = np.array([1,10,1,1])
        elif metrics == 'pr_auc':
            weight = np.array([1,1,10,1])
        elif metrics == 'f1':
            weight = np.array([1,1,1,10])
        elif 'all':
            weight = np.array([1,1,1,1])
        # return negatitive weighted sum (because you want to maximize the sum, 
        # it's equivalent to minimize the negative sum)
        return -np.dot(weight, scores)
  4. Use optimize algorithm differential_evolution() (better then fmin) to find the optimal threshold:

    from scipy import optimize
    output_class_num = Y_test.shape[1]
    bounds = optimize.Bounds([1e-5]*output_class_num,[1]*output_class_num)
    pmodel = proxyModel(model)
    result = optimize.differential_evolution(weighted_score_with_threshold, bounds, args=(pmodel, X_test, Y_test, 'accuracy'))
    # calculate threshold
    threshold = result.x/np.sum(result.x)
    # print the optimized score
    print(scoreFunc(model, X_test, Y_test, threshold_list=threshold))
  • $\begingroup$ you reference @conngrejo who in turn references this soln? are you two related / know each other/friends or just happy coincidence? :) great soln by the way... $\endgroup$
    – JB-007
    Oct 7, 2021 at 0:24

This was helpful, thanks! But it is not applicable during model training. When it comes to using this method after training the model (after finding the hyperparameters relevant to the model), it's valid; only that there has to be some way of standardizing this to avoid loss of generality and for it to be applicable on test data.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.