Multivariate Probability Calibration with Isotonic Bernstein Polynomials

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.

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).