--- title: "Multivariate Probability Calibration with Isotonic Bernstein Polynomials" authors: - "Yongqiao Wang" - "Xudong Liu" doi: "10.24963/ijcai.2020/353" url: "https://ijcai.org/proceedings/2020/353" pdf: "https://ijcai.org/proceedings/2020/0353.pdf" venue: "Proceedings of the Twenty-Ninth International Joint Conference on Artificial Intelligence" year: 2020 pages: "2547-2553" track: "Main track" keywords: - "Machine Learning: Ensemble Methods" --- # Multivariate Probability Calibration with Isotonic Bernstein Polynomials Yongqiao Wang, Xudong Liu *Proceedings of the Twenty-Ninth International Joint Conference on Artificial Intelligence*, pages 2547-2553. DOI: [10.24963/ijcai.2020/353](https://doi.org/10.24963/ijcai.2020/353) ## Abstract Multivariate probability calibration is the problem of predicting class membership probabilities from classification scores of multiple classifiers. To achieve better performance, the calibrating function is often required to be coordinate-wise non-decreasing; that is, for every classifier, the higher the score, the higher the probability of the class labeling being positive. To this end, we propose a multivariate regression method based on shape-restricted Bernstein polynomials. This method is universally flexible: it can approximate any continuous calibrating function with any specified error, as the polynomial degree increases to infinite. Moreover, it is universally consistent: the estimated calibrating function converges to any continuous calibrating function, as the training size increases to infinity. Our empirical study shows that the proposed method achieves better calibrating performance than benchmark methods. ## Full text ## 1 Introduction To predict class membership probability for a given sample is a crucial component in various machine learning-based decision making systems [He et al., 2015; Ustun and Rudin, 2019; Schwarz and Heider, 2019]. One approach to obtain such predictions is to acquire a probability function by estimating the joint distribution function of the class label and the sample vector. However, since the sample vector usually includes both categorical and continuous random variables, estimating this joint distribution function is a challenging task. An alternative is probability calibration that consists of two steps as follows. First, it independently trains M classifiers, whose underlying scoring functions map a feature vector to a dimension-M score vector. Second, it trains a probability calibrating function that maps the dimension-M score vector to a membership probability. We see two advantages of this alternative over direct estimation. One, it incorporates well-established classification algorithms of choice to come up with sample scores. The other, it decreases the dimension of the joint distribution to M, assuming M is much smaller than the dimension of the feature vector. ∗Corresponding Author When M = 1, i.e., there is only one classifier, many parametric and non-parametric methods have been proposed for probability calibration. The related work section of Wang et al. [2019] provides a survey of such methods for binary classification. In general, parametric methods suffer from the lack of flexibility and the vulnerability to mis-specification, while non-parametric calibration methods, including histogram binning and isotonic (or monotonic) regression, usually require a large amount of data to perform reasonably well. Recently, researchers have shown that imposing isotonicity and smoothness on the calibrating function can greatly improve out-of-sample prediction performance, when the training size is small [Zadrozny and Elkan, 2002; Jiang et al., 2011; Naeini and Cooper, 2016]. However, to require monotonicity is a rather difficult task, for it is imposed on every point in the domain along every dimension, which contributes infinitely many inequality constraints. To this end, we propose to solve the multivariate probability calibration problem for binary classified samples, by shape-restricted regression with multivariate Bernstein polynomials (BernPolyFusion, for short). In statistics, univariate Bernstein polynomials have found many successful applications in nonlinear regression under shape-restrictions, such as monotonicity and convexity [Wang and Ghosh, 2012]. The key contributions of this work are as follows. Firstly, this non-parametric method has asymptotic universal flexibility; that is, as the polynomial degree increases, the proposed function family can approximate any continuous multivariate calibrating function with any specified error. This is a great advantage over Ozdemir et al. [2017], in which the specified parametric copula family limits flexibility. Secondly, the fitting calibrating function estimated by this method is coordinate-wise isotonic over the entire domain. This is an advantage over Zhong and Kwok [2013]. In this paper, all deterministic variables are in lower-case letters, while random variables are written in capital letters. Scalars are written in normal letters, and vectors are written in boldfaced letters. For example, the i-th element of constant vector s is si, and the i-th element of random vector X is Xi. We write P[A] the probability of event A, and E[ξ] the mathematical expectation of variable ξ. We denote by I{·} the indicator function, and ≜ reads "equal to by definition". ## 2 ProblemFormulation In a binary classification problem, one considers a random vector (X, Y ), where the feature vector X ∈ X ⊂ Rd and the class label Y ∈ {0, 1}. The classification problem aims to predict the class label Y of each x ∈ X . A typical classification model first estimates a scoring function s : X → R, then predicts the class label of x according to whether its score s(x) is above or below a threshold. A good scoring function is expected to have good ranking power: for any x1, x2 ∈ X , if s(x1) > s(x2), P[Y = 1|X = x1] ≥ P[Y = 1|X = x2]. Without loss of generality, this paper assumes s(X ) ⊂ [0, 1]. Or else one can achieve this by any increasing transformation. Instead of class labels, this paper studies how to predict the conditional probability function Then P[Y = 0|X = x] = 1 - g(x). One can obtain g(·) by estimating the joint distribution of (X, Y ), i.e. two conditional distributions of X|Y = 1 and X|Y = 0, from a training data. However, this task is tough in many applications. An alternative is post-processing that relies on classification models. Assume that M classification models have been independently trained, s(x) ≜ (s1(x), · · · , sM(x))⊤. If the calibrating function is the membership probability for x is f (·) should be estimated by calibration models from the training data {(s(Xi), Yi)}ni=1 . Please note that this calibration is necessary for probability prediction, since many classification methods produce scoring functions that have no direct relationship with membership probability [Guo et al., 2017; Ott et al., 2018; Kuleshov et al., 2018; Kuleshov and Ermon, 2017]. For simplifying notation, we let s ≜ s(x), S ≜ s(X), Si ≜ s(Xi), Dn ≜ {(Si, Yi)}ni=1 . ## 3 RelatedWork ## 3.1 CopulaFusion CopulaFusion, proposed by Ozdemir et al. [2017], estimates f(·) by the Bayesian theorem Each S|Y = ℓ has the following copula representation where Cℓ : [0, 1]M → [0, 1] is the copula function. In this representation, each multivariate distribution consists of M marginal distribution P[Sm ≤ sm|Y = ℓ] and one copula function Cℓ(·). To estimate the calibrating function f(·) from the training data Dn, one should estimate two priors P[Y = ℓ], 2M marginal conditional distributions P[Sm ≤ sm|Y = ℓ], and two copulas Cℓ. P[Y = ℓ] can be replaced with its plug-in estimator ∑ni=1 I{Yi=ℓ}/n. Each P[Sm ≤ sm|Y = ℓ] can be obtained independently with kernel density estimation. Each Cℓ can be obtained by maximum likelihood estimation. ## 3.2 MR-MIC MR-MIC proposed by Zhong and Kwok [2013] estimates the membership probabilities {pi}ni=1 for {Xi}ni=1 with the following optimization where ∑ni=1 [Yi - pi]2 is the empirical error that is a plugin estimator of the L2 risk E[f(S) - Y ]2 , and p⊤Ωp is the manifold regularization for smoothness, and λ > 0 is a hyperparameter for the trade-off between the empirical error and the regularization. Eq. (6b) is the isotonic constraint. E is the transitive reduction of the relation {(i, j)|Si > Sj}. However, the requirement of monotonicity only on ordered pairs of sample points cannot guarantee the monotonicity over the entire domain [0, 1]M . It has the possibility of failing to keep monotone, especially when the training size is small. ## 4 Methodology To keep the presentation simple, we limit our discussion to the case M = 2. The extension to M ≥ 3 is straightforward. This paper estimates the calibrating function f(·) with a shape-restricted non-parametric regression where F is the set where f ∈ C[0,1]2 means that f is continuous over [0, 1]2 . The requirement of coordinate-wise non-decrease arises from the ranking power of each scoring function. The requirement of coordinate-wise non-decrease is very strong. Since it is imposed on every point of [0, 1]2 , the calibrating function is continuously constrained and involves an uncountable number of inequality constraints. Generally, this kind of optimization is a semi-infinite program [Reemtsen and R¨uckmann, 1998], because it involves a finite number of decision variables and an infinite number of inequality constraints. Over-restricted models specify simple functional families that are easy to impose shape restrictions, but these models achieve computational convenience at the expense of flexibility. Under-restricted models simplify these shape restrictions with finitely many inequality constraints on all samples or finite grid points, thus fail to guarantee full adherence. This paper solves the probability calibration problem with shape-restricted multivariate Bernstein polynomials. Univariate Bernstein polynomials are widely regarded as the best shape-preserving functions and have been successfully applied to univariate shape-restricted non-parametric regression [Wang and Ghosh, 2012]. Let f(·) be approximated with a degree-(K1, K2) Bernstein polynomial where bkm,Km (·) is a univariate degree-Km basis Bernstein polynomial BK1K2 is a linear combination of (K1+1)(K2+1) bivariate basis Bernstein polynomials. Proof. This proof is based on three properties of univariate basis Bernstein polynomials Since BK1K2 (0, 0) = β00 and BK1K2 (1, 1) = βK1K2 , we have BK1K2 (0, 0) ≥ 0 ⇔ β00 ≥ 0, and BK1K2 (1, 1) ≤ 1 ⇔ βK1K2 ≤ 1. For all s ∈ [0, 1]2 , we have Because ∀s ∈ [0, 1], bk,K(s) ≥ 0, Eq. (11) implies that BK1K2 (s1, s2) is increasing with respect to s1. The sufficiency of Eq. (12) for s2 is similar. The probability calibration problem with two fore-going classification models can be solved by fitting a bivariate degree-(K1, K2) Bernstein polynomial with training samples {(Si1, Si2, Yi)}ni=1 . The corresponding quadratic program is This quadratic program is convex and tractable, which can be solved efficiently by many off-the-shelf optimization software. If the optimal solution for optimization (13) is β∗ , the estimated calibrating function is let FK be the family of shape-restricted degree-K bivariate Bernstein polynomials ## 4.1 Universal Flexibility This subsection proves that ∪∞K=1FK is dense in F with respect to sup-norm, which means that shape-restricted Bernstein polynomials can be used to approximate any continuous coordinate-wise non-decreasing calibrating function with any arbitrary accuracy as the degree K increases to infinite. Theorem 2. The sequence of function families FK is nested in F, i.e. F1 ⊂ F2 ⊂ · · · ⊂ FK ⊂ · · · ⊂ F ⊂ L2[0, 1]2 , where f ∈ L2[0, 1]2 means ∫[0,1]2 |f(s)|2ds < +∞. Proof. (I) FK ⊂ FK+1, ∀K. For any BK ∈ FK, BK(s) = ∑Kk1=0 ∑Kk2=0 βk1k2 bk1,K(s1)bk2,K(s2), we can rewritten it as a degree-(K + 1) Bernstein polynomial where The above statement is obtained by the iterative property of univariate Bernstein polynomials To verify BK ∈ FK+1, we should further prove that ˜ β ∈ R(K+2)×(K+2) satisfies The first requirement can be verified by ˜ β00 = β00, ˜ β(K+1)(K+1) = βKK and Eq. (10). The second requirement can be obtained by Eq. (11) and The proof for the third requirement is similar. (II) Since all basis Bernstein polynomials and their linear combinations belong to L2[0, 1]2 , FK ⊂ L2[0, 1]2 , ∀K. Theorem 3. ∪∞K=1FK is dense in F with respect to supnorm, i.e. for every f ∈ F, we have Proof. According to DeVore and Lorentz [1993, p.10], for any multivariate continuous function f , there are multivariate Bernstein polynomials that converge uniformly to f as the degrees increase to infinite. So that the set of unconstrained bivariate Bernstein polynomials is dense in C[0,1]2 with respect to super-norm. F ⊂ C[0,1]2 implies that, for each f ∈ F, for any approximation error ϵ > 0, there exists a K ∈ N such that the following function satisfies sups∈[0,1]2 |BK(s)-f(s)| < ϵ. Therefore, if we let β satisfies Eq. (10) - (12) by the properties of f. Thus, BK(·) ∈ FK. Combining this with Theorem 2, we obtain Eq. (19). ## 4.2 Consistency The problem studied in this paper is called random-design regression. Each (Si, Yi) is random, and (S, Y ), (S1, Y1), · · · , (Sn, Yn) are independent and identically distributed. To measure the error of a regression estimate, we use the L2 error: ∫s∈[0,1]2 | ˆ fn(s) - f (s)|2d(s). Since the estimate ˆ fn depends on the data Dn, this L2 error is a random variable. The following theorem on the consistency of the estimator relies on Lemma 1 [Gy¨orfi et al., 2006, Lemma 10.1]. Lemma 1. Let Fn = Fn(Dn) be a class of functions ˜ f : Rd → R depending on the data {(Si, Yi)}ni=1 . If ˆ fn is obtained by then where µ denotes the distribution of S. Theorem 4. Provided that ˆ fn is obtained by Eq. (13) - (14), if the degree K satisfies then, (I) ˆ fn is strongly universally consistent, i.e. for any f ∈ F, (II) ˆ fn is weakly universally consistent, i.e. for any f ∈ F, as n → ∞, Proof. (I) STRONG CONSISTENCY. Due to the Lemma 1, to prove the strong universal consistency of ˆ fn, it suffices to verify that, as n →∞, Let us first prove Eq. (27). According to Theorem 3, ∪∞K=1FK is dense in F with respect to sup norm, which follows that ∪∞K=1FK is dense in F with respect to L2(µ) for any distribution µ. It implies that, for every f ∈ F, for any arbitrarily small error ϵ > 0, there exists a BK0(ϵ) ∈ ∪∞K=1FK that satisfies ∫[0,1]2 |BK0(ϵ)(s) - f(s)|2µ(ds) < ϵ. Because K ↑ ∞ as n → ∞, there is a n0(ϵ) ∈ N such that K > K0(ϵ) for all n ≥ n0(ϵ). Combining this with the nested property of Fns, one concludes that Since ϵ > 0 is arbitrary, this implies Eq. (27). To verify Eq. (26), let F + n be the class of all subgraphs of functions of Fn and VF + n be the Vapnik-Chervonenkis (VC) dimension of F + n . Since each ˜ f ∈ Fn is a linear combination of (K + 1)2 basis Bernstein polynomials, and F + n is a subset of which is a vector space with dimension (K + 1)2 + 1 of real functions on R3 , by Theorem 9.5 of [Gy¨orfi et al., 2006], we have The proof for the first inequality is same as [Wang et al., 2019, Theorem 4]. The second inequality follows from log(x) ≤ x - 1 ≤ x for x > 0. Therefore, Hence, provided that limn→∞K2 /n → 0, we have By the Borel-Cantelli lemma, we obtain Eq. (26). Therefore, Eq. (23) is a sufficient condition for the strong universal convergence of ˆ fn(·). (II) WEAK CONSISTENCY. To prove the weak universal consistency of ˆ fn(·), it suffices to verify, as n→∞, Because Eq. (34) directly follows fromEq. (28), it is sufficient to verify Eq. (33). If ξ is a nonnegative random variable, for an arbitrary ϵ > 0, we have Following this, The first inequality is obtained by Eq. (35) and (31). Eq. (33) holds if K2 /n → 0 as n →∞. Hence, the condition (23) is sufficient for the weak universal consistency of ˆ fn(·). ## 5 Experiments In this experiment, to make each sm(X ) ∈ [0, 1], every score sm(Xi) is transformed by the empirical cumulative distribution function using {sm(Xi)}ni=1 . In CopulaFusion, the copula function for each class is chosen according to the criterion of maximum likelihood. Candidate copula families are Gaussian, Clayton, rotated Clayton, Plackett, Frank, Gumbel, rotated Gumbel and Student's t. In MR-MIC, Ω is the similarity matrix with ωij = 1/∥Si -Sj∥. Because model (Eq.6) only obtains calibrated probabilities {ˆ pi}ni=1 for training samples {Si}ni=1 , we estimate the calibrating function by linear interpolation with the scattered data {(Si, ˆ pi)}ni=1 . Quadratic optimizations for MR-MIC and BernPolyFusion are solved with CVX [Grant and Boyd, 2014]. ## 5.1 Characteristics of Calibrating Functions The data is SUSY from UCI Machine Learning Repository [Dua and Graff, 2017]. Two fore-going classifiers are feedforward networks with hidden layer sizes 10 and 20. We randomly draw 200 samples for training two classifiers and 400 samples for training the calibrating model. The hyperparameters are arbitrarily chosen: λ = 10-5 in MR-MIC, and K1 = K2 = 5 in the proposed method. Figure 1(a) shows that both classifiers have good ranking power. Figure 1(b)-1(d) shows the main characteristics of three estimated calibrating functions with contours. In all figures, by and large, the estimated calibrating probability increases as S moves from southwest to northeast. Among the three methods, only BernPolyFusion can adhere to the requirement of coordinate-wise non-decrease. MR-MIC imposes this requirement only on training samples {Si}ni=1 , while CopulaFusion does not explicitly impose this requirement. The calibrating function from BernPolyFusion is more smooth than MR-MIC and CopulaFusion. ## 5.2 Model Comparison Ideal performance measures for probability calibration should be based on the difference between the true f and the estimated ˆ f, e.g. sups∈[0,1] |f(s)- ˆ f(s)| and ∫s∈[0,1] |f(s)- ˆ f(s)|ds. However, in practice the true f is unavailable. When the size of the test data is very large, f can be approximated in the following way: the domain [0, 1]2 is discretized with a uniformly-spaced grid (10×10) and the true probability at each grid cell is approximated with the percent of class-1 among samples that falls in this grid cell. Because the majority of scores scatter around the antidiagonal, there are few test samples scatter in grid cells around northwest and southeast corners. Thus, approximated probabilities at these grid cells are subject to large deviations. Therefore, each grid cell with the number of test samples smaller than 1000 is discarded from performance measurement. The following two measures are used: MCE = maxs∈G |f(s) - ˆ f(s)| and ECE = ∑s∈G |f(s) - ˆ f (s)|/|G|, where G is the set of centers of grid cells with the test number no less than 1000. The third performance measure is the Brier score, which is the mean squared difference between predicted probabilities and class labels. Figure 1: Training data and contours of three calibrating functions ![Image](https://ijcai.org/proceedings/2020/fulltext/353/image_000000_00dcaae24c8d920a8831b1e3dc1f0e22fd6ce2262253cd1fa3df15d1eb6f1913.png) Table 1: Model Comparison (%. The best performance is highlighted.) | | Method | Adult | Census | Covertype | Dota2 | SUSY | Mean | |-------|----------------|--------------|--------------|--------------|--------------|--------------|--------------| | MCE | MR-MIC | 6.142±0.843 | 7.116±0.825 | 5.785±0.611 | 5.353±0.576 | 6.185±0.479 | 6.116±0.667 | | MCE | CopulaFusion | 5.309±0.723 | 6.817±0.611 | 5.882±0.626 | 6.290±0.508 | 6.592±0.469 | 6.178±0.587 | | MCE | BernPolyFusion | 5.613±0.650 | 6.241±0.565 | 5.367±0.434 | 5.591±0.528 | 5.054±0.394 | 5.573±0.514 | | ECE | MR-MIC | 2.459±0.441 | 3.787±0.468 | 3.444±0.359 | 3.124±0.393 | 3.182±0.258 | 3.199±0.384 | | ECE | CopulaFusion | 2.548±0.296 | 3.890±0.444 | 3.369±0.326 | 3.319±0.343 | 2.963±0.240 | 3.218±0.330 | | ECE | BernPolyFusion | 2.226±0.309 | 3.587±0.370 | 2.847±0.231 | 2.799±0.268 | 2.432±0.241 | 2.727±0.284 | | Brier | MR-MIC | 15.485±0.915 | 19.719±1.338 | 19.776±0.861 | 18.740±0.758 | 19.015±0.646 | 18.547±0.904 | | Brier | CopulaFusion | 16.441±0.703 | 19.714±1.095 | 18.067±0.824 | 19.264±0.711 | 17.771±0.506 | 18.252±0.768 | | Brier | BernPolyFusion | 14.375±0.854 | 18.966±0.724 | 16.095±0.577 | 16.811±0.434 | 15.539±0.473 | 16.357±0.612 | Experiments are performed on five large data from [Dua and Graff, 2017] with size larger than 40,000: Adult, Census, Covertype, Dota2 and SUSY. In Covertype, the largest class (Lodgepole Pine) is the positive class and others are the negative class. For each data, three size-500 datasets are randomly drawn. The first dataset is used for training classifier 1, and the second dataset is used for training classifier 2, and the third dataset is used for training calibration models. All other samples are used for testing calibration models. In each training step, we use 5-fold cross-validation to determine hyperparameters. Two fore-going classifiers are three-layer feedforward networks. The above procedure is repeated 10 times, and model comparison is based on the average and standard deviation of ten rounds. Table 1 demonstrates that the proposed method outperforms benchmark models in all but two instances, where CopulaFusion bests on Adult and MR-MIC on Dota2, both with respect to MCE measure. Considering the results averaged over all five datasets, we see our BernPolyFusion unanimously performs the best. ## 5.3 Running Time The complexity of CopulaFusion depends on the choice of the copula family. For some copula families, the optimization for the maximum likelihood estimation is non-convex. The training of MR-MIC involves a quadratic program with n decision variables. An ADMM-based algorithm for this optimization has a complexity O(n2 ). The training of BernPolyFusion involves a quadratic program with (K + 1)2 decision variables and 2K(K + 1) + 2 constraints. Both numbers are independent of the training size n. According to Section 4.6.2 of [Ben-Tal, 2019], Data(P) and Size(p) are linearly increasing with n, thus the Newton complexity of ϵ-solution is O(n). For the aforementioned experiments in Subsection 5.2 where n = 500, the CPU time in seconds for training MRMIC, CopulaFusion, and BernPolyFusion are 1.335 ± 0.128, 1.486±0.172, and 0.501±0.044, respectively. If n is 1,000, the CPU time in seconds for three methods are 7.964±0.807, 1.756 ± 0.246 and 0.872 ± 0.090, respectively. Therefore, compared to MR-MIC and CopulaFusion, we see that BernPolyFusion attains less running time and better scalability. ## 6 Conclusions A novel method based on shape-restricted Bernstein polynomial regression is proposed for probability calibration. This method has universal fitting flexibility and is both strongly and weakly universally consistent. Experimental results show that this method has a great advantage over benchmark methods. Future work includes the extension from binary-class to multi-class or structured classification problems [Kuleshov and Liang, 2015; Leathart et al., 2019]. ## Acknowledgments The work of the first author was supported by NSFC (71571163) and Zhejiang NSF (LY19G010001). ## References 1. [Ben-Tal, 2019] Aharon Ben-Tal. Lectures on modern con- vex optimization. Technical report, School of Industrial & Systems Engineering, Georgia Institute of Technology, 2019. 2. [DeVore and Lorentz, 1993] Ronald A. DeVore and George G. Lorentz. Constructive Approximation. Springer-Verlag, 1993. 3. [Dua and Graff, 2017] Dheeru Dua and Casey Graff. UCI Machine Learning Repository, 2017. 4. [Grant and Boyd, 2014] Michael Grant and Stephen Boyd. CVX: Matlab software for disciplined convex program- ming, version 2.1. http://cvxr.com/cvx, March 2014. 5. [Guo et al., 2017] Chuan Guo, Geoff Pleiss, Yu Sun, and K- ilian Q. Weinberger. On calibration of modern neural net- works. In Proc. of ICML'17, pages 1321-1330, 2017. 6. [Györfi et al., 2006] László Györfi, Michael Kohler, Adam Krzyzak, and Harro Walk. A Distribution-free Theory of Nonparametric Regression. Springer Science & Business Media, 2006. 7. [He et al., 2015] Jiazhen He, James Bailey, Benjamin I. P. Rubinstein, and Rui Zhang. Identifying at-risk students in massive open online courses. In Proc. of AAAI'15, pages 1749-1755, 2015. 8. [Jiang et al., 2011] Xiaoqian Jiang, Melanie Osl, Jihoon Kim, and Lucila Ohno-Machado. Smooth isotonic regres- sion: A new method to calibrate predictive models. AMIA Summits on Translational Science Proceedings, 2011:16, 2011. 9. [Kuleshov and Ermon, 2017] Volodymyr Kuleshov and Ste- fano Ermon. Estimating uncertainty online against an ad- versary. In Proc. of AAAI'17, pages 2110-2116, 2017. 10. [Kuleshov and Liang, 2015] Volodymyr Kuleshov and Per- cy S Liang. Calibrated structured prediction. In Ad- vances in Neural Information Processing Systems 28, pages 3474-3482. MIT Press, 2015. 11. [Kuleshov et al., 2018] Volodymyr Kuleshov, Nathan Fen- ner, and Stefano Ermon. Accurate uncertainties for deep learning using calibrated regression. In Proc. of ICML'18, pages 2801-2809, 2018. 12. [Leathart et al., 2019] Tim Leathart, Eibe Frank, Bernhard Pfahringer, and Geoffrey Holmes. On calibration of nested dichotomies. In Proc. of PAKDD'19, pages 69-80, 2019. 13. [Naeini and Cooper, 2016] Mahdi Pakdaman Naeini and Gregory F Cooper. Binary classifier calibration using an ensemble of near isotonic regression models. In Proc. of ICDM'16, pages 360-369, 2016. 14. [Ott et al., 2018] Myle Ott, Michael Auli, David Grangier, and Marc'Aurelio Ranzato. Analyzing uncertainty in neu- ral machine translation. In Proc. of ICML'18, pages 3956- 3965, 2018. 15. [Ozdemir et al., 2017] Onur Ozdemir, Thomas G Allen, So- ra Choi, Thakshila Wimalajeewa, and Pramod K Varsh- ney. Copula based classifier fusion under statistical depen- dence. IEEE Transactions on Pattern Analysis and Ma- chine Intelligence, 40(11):2740-2748, 2017. 16. [Reemtsen and Rückmann, 1998] Rembert Reemtsen and Jan-J Rückmann. Semi-infinite Programming, volume 25. Springer Science & Business Media, 1998. 17. [Schwarz and Heider, 2019] Johanna Schwarz and Dominik Heider. GUESS: projecting machine learning scores to well-calibrated probability estimates for clinical decision- making. Bioinformatics, 35(14):2458-2465, 2019. 18. [Ustun and Rudin, 2019] Berk Ustun and Cynthia Rudin. Learning optimized risk scores. Journal of Machine Learning Research, 20:1-75, 2019. 19. [Vellanki et al., 2019] Pratibha Vellanki, Santu Rana, Sunil Gupta, David Rubin de Celis Leal, Alessandra Sutti, Mur- ray Height, and Svetha Venkatesh. Bayesian functional optimisation with shape prior. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 33, pages 1617-1624, 2019. 20. [Wang and Ghosh, 2012] J. Wang and S.K. Ghosh. Shape restricted nonparametric regression with bernstein poly- nomials. Computational Statistics & Data Analysis, 56(9):2729 -2741, 2012. 21. [Wang et al., 2019] Y. Wang, L. Li, and C. Dang. Calibrat- ing classification probabilities with shape-restricted poly- nomial regression. IEEE Transactions on Pattern Analysis and Machine Intelligence, 41(8):1813-1827, 2019. 22. [Zadrozny and Elkan, 2002] Bianca Zadrozny and Charles Elkan. Transforming classifier scores into accurate mul- ticlass probability estimates. In Proc. of KDD'02, pages 694-699, 2002. 23. [Zhong and Kwok, 2013] Leon Wenliang Zhong and James T Kwok. Accurate probability calibration for mul- tiple classifiers. In Proc. of IJCAI'13, pages 1939-1945, 2013.