import math
import numpy as np
import scipy

class BaseConformalPredictor:
    def __init__(self, label_type, nonconformity_score_type):
        # The label and non-conformity score types (continuous or discrete)
        self.label_type = label_type
        self.nonconformity_score_type = nonconformity_score_type

    def set_train(self, X_train, y_train):
        # Train the machine learning model
        self.model.fit(X_train, y_train)

    def set_calibration(self, X_cal, y_cal, alpha):
        # Compute the variables dependent on the calibration dataset size and
        # the significance level
        self.n = X_cal.shape[0]
        self.n_alpha = math.ceil((1. - alpha) * (self.n + 1)) - 1

        # Compute the conformal threshold
        r_cal = self.get_nonconformity_score(X_cal, y_cal)
        r_cal.sort()
        r = np.append(r_cal, np.inf)
        self.threshold = r[self.n_alpha]

    def get_prediction_set_size_estimates_known_multiplicative_factor(
        self, X_access, y_access, gamma
    ):
        # Compute the variables dependent on the accessible dataset size
        k = X_access.shape[0]
        delta = math.sqrt(math.log(2. / gamma) / (2. * k))

        # Compute the accessible dataset non-conformity scores
        r_access = self.get_nonconformity_score(X_access, y_access)

        # The space of non-conformity scores over which to integrate/sum
        if self.nonconformity_score_type == 'continuous':
            # The accessible dataset non-confomrity scores, with r_min and
            # r_max, because the empirical distribution of non-conformity scores
            # is a step function over this space
            r_access.sort()
            r = np.insert(r_access, 0, self.r_min)
            r = np.append(r, self.r_max)
        elif self.nonconformity_score_type == 'discrete':
            # The space of non-conformity scores
            r = self.r_space

        # Compute the desired epirical distribution of non-conformity scores
        p = (r_access[np.newaxis, :] < r[:, np.newaxis]).mean(1)
        p_lower = (p - delta).clip(0., 1.)
        p_upper = (p + delta).clip(0., 1.)

        # Compute the binomial CDF terms
        binom_cdfs = scipy.stats.binom.cdf(
            self.n_alpha, self.n, [p, p_upper, p_lower]
        )

        # Compute the expected prediction set size estimates
        if self.nonconformity_score_type == 'continuous':
            # Compute the integral over the non-conformity score space
            int_mf = self.integrate_multiplicative_factor(r)
            sizes = (binom_cdfs[:, 1 :] * int_mf[np.newaxis, :]).sum(1)
        elif self.nonconformity_score_type == 'discrete':
            # Compute the sum over the non-conformity score space
            mf = self.multiplicative_factor(r)
            sizes = (binom_cdfs * mf[np.newaxis, :]).sum(1)

        # Return the expected prediction set size estimates
        return sizes

    def get_prediction_set_size_estimates_unknown_multiplicative_factor(
        self, X_access, y_access, gamma
    ):
        # Compute the variables dependent on the accessible dataset size
        k = X_access.shape[0]
        delta = math.sqrt(math.log(2. / gamma) / (2. * k))

        # Compute the accessible dataset non-conformity scores
        r_access = self.get_nonconformity_score(X_access, y_access)

        # Compute the accessible dataset non-conformity scores for all labels
        R_access = self.get_nonconformity_scores(X_access)

        # Compute the desired epirical distribution of non-conformity scores
        p = (
            r_access[np.newaxis, np.newaxis, :] < R_access[:, :, np.newaxis]
        ).mean(2)
        p_lower = (p - delta).clip(0., 1.)
        p_upper = (p + delta).clip(0., 1.)

        # Compute the binomial CDF terms
        binom_cdfs = scipy.stats.binom.cdf(
            self.n_alpha, self.n, [p, p_upper, p_lower]
        )

        # Compute the expected prediction set size estimates
        if self.label_type == 'continuous':
            # Compute the mean over the accessible dataset, and the integral
            # over the label space
            sizes = self.y_width * binom_cdfs.mean(1).sum(1)
        elif self.label_type == 'discrete':
            # Compute the mean over the accessible dataset, and the sum over the
            # label space
            sizes = binom_cdfs.mean(1).sum(1)

        # Return the expected prediction set size estimates
        return sizes


base_conformal_predictor.py
Guneet-Dhillon/expected-conformal-prediction-set-size
import numpy as np

from sklearn.ensemble import RandomForestClassifier, RandomForestRegressor
from quantile_forest import RandomForestQuantileRegressor

from base_conformal_predictor import BaseConformalPredictor

class LPRegressionConformalPredictor(BaseConformalPredictor):
    def __init__(self, y, p=1., model='randomforest'):
        super().__init__('continuous', 'continuous')

        # Set the order of the norm
        self.p = p

        # Set the minimum and the maximum non-conformity scores
        self.r_min = 0.
        self.r_max = (y.max() - y.min()) ** p

        # The machine learning model
        if model == 'randomforest':
            self.model = RandomForestRegressor()

        # The multiplicative factor is known
        self.get_prediction_set_size_estimates = \
            self.get_prediction_set_size_estimates_known_multiplicative_factor

    def get_model_prediction(self, X):
        # Compute the machine learning model prediction
        yh = self.model.predict(X)
        return yh

    def get_nonconformity_score(self, X, y):
        # Compute the non-conformity score
        yh = self.get_model_prediction(X)
        r = np.abs(yh - y) ** self.p
        return r

    def get_prediction(self, X_test):
        # Compute the conformal prediction set
        yh_test = self.get_model_prediction(X_test)
        C = [
            yh_test - (self.threshold ** (1. / self.p)),
            yh_test + (self.threshold ** (1. / self.p))
        ]
        return C

    def integrate_multiplicative_factor(self, r):
        # Compute the integral of the multiplicative factor over each
        # (r[i], r[i + 1]] region
        int_mf = 2. * ((r[1 :] ** (1. / self.p)) - (r[: -1] ** (1. / self.p)))
        return int_mf

class ZeroOneClassificationConformalPredictor(BaseConformalPredictor):
    def __init__(self, num_classes, model='randomforest'):
        super().__init__('discrete', 'discrete')

        # Set the number of classes
        self.num_classes = num_classes

        # Set the space of non-conformity scores
        self.r_space = np.array([0., 1.])

        # The machine learning model
        if model == 'randomforest':
            self.model = RandomForestClassifier()

        # The multiplicative factor is known
        self.get_prediction_set_size_estimates = \
            self.get_prediction_set_size_estimates_known_multiplicative_factor

    def get_model_prediction(self, X):
        # Compute the machine learning model prediction
        yh = self.model.predict(X)
        return yh

    def get_nonconformity_score(self, X, y):
        # Compute the non-conformity score
        yh = self.get_model_prediction(X)
        r = (yh != y).astype(float)
        return r

    def get_prediction(self, X_test):
        # Compute the conformal prediction set
        n_test = X_test.shape[0]
        if self.threshold == 0.:
            yh_test = self.get_model_prediction(X_test)
            C = np.zeros((n_test, self.num_classes), bool)
            C[np.arange(n_test), yh_test] = 1.
        else:
            C = np.ones((n_test, self.num_classes), bool)
        return C

    def multiplicative_factor(self, r):
        # Compute the multiplicative factor for each non-conformity score
        mf = (r == 0.) + ((self.num_classes - 1) * (r == 1.))
        return mf

class CQRRegressionConformalPredictor(BaseConformalPredictor):
    def __init__(self, y, y_num=100, model='randomforest'):
        super().__init__('continuous', 'continuous')

        # Set the minimum and the maximum labels
        self.y_min = y.min()
        self.y_max = y.max()

        # Set the number of labels to use for the expected prediction set size
        # estimates
        self.y_num = y_num

        # The machine learning model
        if model == 'randomforest':
            self.model = RandomForestQuantileRegressor()

        # The multiplicative factor is unknown
        self.get_prediction_set_size_estimates = \
            self.get_prediction_set_size_estimates_unknown_multiplicative_factor

    def set_calibration(self, X_cal, y_cal, alpha):
        # Set the significance level
        self.alpha = alpha
        super().set_calibration(X_cal, y_cal, alpha)

    def get_model_prediction(self, X):
        # Compute the machine learning model prediction
        quantiles = [self.alpha / 2., 1. - (self.alpha / 2.)]
        Yh = self.model.predict(X, quantiles=quantiles)
        return Yh

    def get_nonconformity_score(self, X, y):
        # Compute the non-conformity score
        Yh = self.get_model_prediction(X)
        r = np.maximum(Yh[:, 0] - y, y - Yh[:, 1])
        return r

    def get_nonconformity_scores(self, X):
        # Compute all labels to use
        ys, self.y_width = np.linspace(
            self.y_min, self.y_max, self.y_num + 1, retstep=True
        )
        ys += (self.y_width / 2.)
        ys = ys[: -1]

        # Compute the non-conformity scores for all labels
        Yh = self.get_model_prediction(X)
        R = []
        for y in ys:
            r = np.maximum(Yh[:, 0] - y, y - Yh[:, 1])
            R.append(r[:, np.newaxis])
        R = np.hstack(R)
        return R

    def get_prediction(self, X_test):
        # Compute the conformal prediction set
        Yh_test = self.get_model_prediction(X_test)
        C = [Yh_test[:, 0] - self.threshold, Yh_test[:, 1] + self.threshold]
        return C

class LACClassificationConformalPredictor(BaseConformalPredictor):
    def __init__(self, num_classes, model='randomforest'):
        super().__init__('discrete', 'continuous')

        # Set the number of classes
        self.num_classes = num_classes

        # The machine learning model
        if model == 'randomforest':
            self.model = RandomForestClassifier()

        # The multiplicative factor is unknown
        self.get_prediction_set_size_estimates = \
            self.get_prediction_set_size_estimates_unknown_multiplicative_factor

    def get_model_prediction(self, X):
        # Compute the machine learning model prediction
        Yh = self.model.predict_proba(X)
        return Yh

    def get_nonconformity_score(self, X, y):
        # Compute the non-conformity score
        n = X.shape[0]
        Yh = self.get_model_prediction(X)
        r = 1. - Yh[np.arange(n), y]
        return r

    def get_nonconformity_scores(self, X):
        # Compute the non-conformity scores for all labels
        Yh = self.get_model_prediction(X)
        R = 1. - Yh
        return R

    def get_prediction(self, X_test):
        # Compute the conformal prediction set
        R_test = self.get_nonconformity_scores(X_test)
        C = R_test <= self.threshold
        return C

class APSClassificationConformalPredictor(BaseConformalPredictor):
    def __init__(self, num_classes, model='randomforest'):
        super().__init__('discrete', 'continuous')

        # Set the number of classes
        self.num_classes = num_classes

        # The machine learning model
        if model == 'randomforest':
            self.model = RandomForestClassifier()

        # The multiplicative factor is unknown
        self.get_prediction_set_size_estimates = \
            self.get_prediction_set_size_estimates_unknown_multiplicative_factor

    def get_model_prediction(self, X):
        # Compute the machine learning model prediction
        Yh = self.model.predict_proba(X)
        return Yh

    def get_nonconformity_score(self, X, y):
        # Compute the non-conformity score
        n = X.shape[0]
        Yh = self.get_model_prediction(X)
        yh = Yh[np.arange(n), y]
        r = \
            ((yh[:, np.newaxis] < Yh) * Yh).sum(1) \
            + (np.random.uniform(size=n) * yh)
        return r

    def get_nonconformity_scores(self, X):
        # Compute the non-conformity scores for all labels
        n = X.shape[0]
        Yh = self.get_model_prediction(X)
        R = []
        for y in range(self.num_classes):
            yh = Yh[:, y]
            r = \
                ((yh[:, np.newaxis] < Yh) * Yh).sum(1) \
                + (np.random.uniform(size=n) * yh)
            R.append(r[:, np.newaxis])
        R = np.hstack(R)
        return R

    def get_prediction(self, X_test):
        # Compute the conformal prediction set
        R_test = self.get_nonconformity_scores(X_test)
        C = R_test <= self.threshold
        return C


conformal_predictors.py
Guneet-Dhillon/expected-conformal-prediction-set-size
import numpy as np
import os
import pickle

from tqdm import tqdm

from utils import get_size_error, get_mean_std

def get_stats():
    stats = {}
    for file_name in tqdm(sorted(os.listdir('./results')), leave=False):
        # Get the saved results
        file_path = os.path.join('./results', file_name)
        file = open(file_path, 'rb')
        stat_in = pickle.load(file)
        file.close()

        # Get the saved stats
        typ = stat_in['type']
        error = np.array(stat_in['results']['error'])
        size = np.array(stat_in['results']['size'])
        size_est = np.array(stat_in['results']['estimated size'])
        size_est_lower = np.array(stat_in['results']['estimated size (lower)'])
        size_est_upper = np.array(stat_in['results']['estimated size (upper)'])

        # Compute the mean and the standard deviation of the stats
        error = get_mean_std(error)
        size = get_mean_std(size)

        # Compute the prediction set size interval error
        error_size_est = get_size_error(size[0], size_est_lower, size_est_upper)

        # Compute the mean and the standard deviation of the stats
        size_est = get_mean_std(size_est)
        size_est_lower = get_mean_std(size_est_lower)
        size_est_upper = get_mean_std(size_est_upper)
        error_size_est = get_mean_std(error_size_est)

        # Accumulate the results
        stat_out = {
            'dataset'                           : stat_in['dataset'],
            'num_datapoints'                    : stat_in['num_datapoints'],
            'num_features'                      : stat_in['num_features'],
            'error'                             : error,
            'size'                              : size,
            'estimated size'                    : size_est,
            'estimated size (lower)'            : size_est_lower,
            'estimated size (upper)'            : size_est_upper,
            'error (estimated size interval)'   : error_size_est
        }
        if 'Regression' in typ:
            stat_out['y_range'] = stat_in['y_range']
        elif 'Classification' in typ:
            stat_out['num_classes'] = stat_in['num_classes']

        # Sort the results by the conformal predictor type
        if typ not in stats:
            stats[typ] = []
        stats[typ].append(stat_out)

    return stats

def print_stats(stats):
    for typ in stats.keys():
        print('Type : %s' % typ)
        print()

        # Print the dataset stats
        print('Dataset')
        if 'Regression' in typ:
            print(
                'Dataset \t Number of data points \t Number of features '
                '\t Label range'
            )
            for stat in stats[typ]:
                print('%s \t %d \t %d \t [%.4f, %.4f]' % (
                    stat['dataset'], stat['num_datapoints'],
                    stat['num_features'], *stat['y_range']
                ))
        elif 'Classification' in typ:
            print(
                'Dataset \t Number of data points \t Number of features '
                '\t Number of classes'
            )
            for stat in stats[typ]:
                print('%s \t %d \t %d \t %d' % (
                    stat['dataset'], stat['num_datapoints'],
                    stat['num_features'], stat['num_classes']
                ))
        print()

        # Print the conformal prediction coverage stats
        print('Coverage')
        print('Dataset \t Error')
        for stat in stats[typ]:
            print('%s \t %.4f' % (stat['dataset'], stat['error'][0]))
        print()

        # Print the conformal prediction set size stats
        print('Prediction set size')
        print(
            'Dataset \t Interval lower bound \t Observed \t Estimated '
            '\t Interval upper bound \t Interval error'
        )
        for stat in stats[typ]:
            print(
                '%s \t %.2f (%.2f) \t %.2f (%.2f) \t %.2f (%.2f) '
                '\t %.2f (%.2f) \t %.4f' % (
                    stat['dataset'], *stat['estimated size (lower)'],
                    *stat['size'], *stat['estimated size'],
                    *stat['estimated size (upper)'],
                    stat['error (estimated size interval)'][0]
            ))
        print()

def main():
    stats = get_stats()
    print_stats(stats)

if __name__ == '__main__':
    main()


main_analyze.py
Guneet-Dhillon/expected-conformal-prediction-set-size
import argparse
import json
import numpy as np
import os
import pickle

from tqdm import tqdm

from conformal_predictors import \
    LPRegressionConformalPredictor, ZeroOneClassificationConformalPredictor, \
    CQRRegressionConformalPredictor, LACClassificationConformalPredictor, \
    APSClassificationConformalPredictor
from utils import set_random_seed, get_datasets, get_error, get_size

def get_dataset_results(dataset, X, y, args):
    typ = args['type']
    alpha = args['alpha']
    gamma = args['gamma']
    frac_train = args['frac_train']
    frac_cal = args['frac_cal']
    it_train = args['it_train']
    it_cal = args['it_cal']

    # Set random seeds
    set_random_seed()

    # Compute dataset sizes
    n = X.shape[0]
    n_train = round(frac_train * n)
    n_cal = round(frac_cal * n)
    n_test = n - n_train - n_cal

    # Initialize the conformal predictor
    if typ == 'LPRegression':
        predictor = LPRegressionConformalPredictor(y)
    elif typ == 'ZeroOneClassification':
        num_classes = np.unique(y).shape[0]
        predictor = ZeroOneClassificationConformalPredictor(num_classes)
    elif typ == 'CQRRegression':
        predictor = CQRRegressionConformalPredictor(y)
    elif typ == 'LACClassification':
        num_classes = np.unique(y).shape[0]
        predictor = LACClassificationConformalPredictor(num_classes)
    elif typ == 'APSClassification':
        num_classes = np.unique(y).shape[0]
        predictor = APSClassificationConformalPredictor(num_classes)

    results = {
        'error' : [], 'size' : [], 'estimated size' : [],
        'estimated size (lower)' : [], 'estimated size (upper)' : []
    }

    for _ in tqdm(range(it_train), leave=False):
        idx = np.arange(n)

        # Check the presence of all classes in the classification case
        while True:
            # Randomly sample the training dataset
            np.random.shuffle(idx)
            idx_train = idx[: n_train]
            idx_no_train = idx[n_train :]

            X_train = X[idx_train]
            y_train = y[idx_train]

            if 'Regression' in typ:
                break
            elif 'Classification' in typ:
                if num_classes == np.unique(y_train).shape[0]:
                    break

        # Set the training dataset
        predictor.set_train(X_train, y_train)

        for _ in tqdm(range(it_cal), leave=False):
            # Randomly sample the calibration dataset
            np.random.shuffle(idx_no_train)
            idx_cal = idx_no_train[: n_cal]
            idx_no_train_cal = idx_no_train[n_cal :]

            X_cal = X[idx_cal]
            y_cal = y[idx_cal]

            # Set the calibration dataset
            predictor.set_calibration(X_cal, y_cal, alpha)

            # Compute the estimates for the expected prediction set size
            X_access = X_cal
            y_access = y_cal
            size, size_lower, size_upper = \
                predictor.get_prediction_set_size_estimates(
                    X_access, y_access, gamma
                )
            results['estimated size'].append(size)
            results['estimated size (lower)'].append(size_lower)
            results['estimated size (upper)'].append(size_upper)

            # Sample the remaining test dataset
            idx_test = idx_no_train_cal

            X_test = X[idx_test]
            y_test = y[idx_test]

            # Compute the prediction set and its corresponding error and size
            C = predictor.get_prediction(X_test)
            error = get_error(typ, C, y_test)
            size = get_size(typ, C)
            results['error'].append(error.tolist())
            results['size'].append(size.tolist())

        # Accumulate the results
        stat = {
            'type'              : typ,
            'dataset'           : dataset,
            'alpha'             : alpha,
            'gamma'             : gamma,
            'num_datapoints'    : n,
            'num_features'      : X.shape[1],
            'num_train'         : n_train,
            'num_cal'           : n_cal,
            'num_test'          : n_test,
            'it_train'          : it_train,
            'it_cal'            : it_cal,
            'results'           : results
        }
        if 'Regression' in typ:
            stat['y_range'] = [y.min(), y.max()]
        elif 'Classification' in typ:
            stat['num_classes'] = num_classes

        # Save the results
        os.makedirs('results', exist_ok=True)
        name = {'type' : typ, 'dataset' : dataset}
        file_path = os.path.join('./results', '%s.pkl' % json.dumps(name))
        file = open(file_path, 'wb')
        pickle.dump(stat, file)
        file.close()

def get_results(args):
    # Run each dataset individually
    for dataset, X, y in tqdm(get_datasets(args['type']), leave=False):
        get_dataset_results(dataset, X, y, args)

def main():
    # Get inputs
    parser = argparse.ArgumentParser(
        description='On the Expected Size of Conformal Prediction Sets'
    )
    parser.add_argument(
        '--type', type=str, required=True,
        choices=[
            'LPRegression', 'ZeroOneClassification', 'CQRRegression',
            'LACClassification', 'APSClassification'
        ],
        help='Conformal predictor type'
    )
    parser.add_argument(
        '--alpha', type=float, default=0.1,
        help='Conformal predictor significance level'
    )
    parser.add_argument(
        '--gamma', type=float, default=0.1,
        help='Prediction set size interval significance level'
    )
    parser.add_argument(
        '--frac_train', type=float, default=0.25,
        help='Fraction used as training dataset'
    )
    parser.add_argument(
        '--frac_cal', type=float, default=0.25,
        help='Fraction used as calibration dataset'
    )
    parser.add_argument(
        '--it_train', type=int, default=10,
        help='Number of iterations for sampling a new training dataset'
    )
    parser.add_argument(
        '--it_cal', type=int, default=100,
        help='Number of iterations for sampling a new calibration dataset'
    )
    args = vars(parser.parse_args())

    # Run the main code
    get_results(args)

if __name__ == '__main__':
    main()


main_run.py
Guneet-Dhillon/expected-conformal-prediction-set-size
import math
import numpy as np
import random

from ucimlr import classification_datasets, regression_datasets

def set_random_seed(seed=42):
    # Set random seeds
    random.seed(seed)
    np.random.seed(seed)

def get_regression_datasets():
    # Get UCI regression datasets
    datasets = []
    for dataset in regression_datasets.all_datasets():
        train = dataset('data', split='train', validation_size=0.)
        test = dataset('data', split='test', validation_size=0.)

        name = train.name
        X = np.vstack([train.x, test.x])
        y = np.vstack([train.y, test.y])

        if X.shape[0] >= 1000 and y.shape[1] == 1:
            y = y.squeeze()

            dataset = (name, X, y)
            datasets.append(dataset)

    return datasets

def get_classification_datasets():
    # Get UCI classification datasets
    datasets = []
    for dataset in classification_datasets.all_datasets():
        train = dataset('data', split='train', validation_size=0.)
        test = dataset('data', split='test', validation_size=0.)

        name = train.name
        X = np.vstack([train.x, test.x])
        y = np.hstack([train.y, test.y])

        if X.shape[0] >= 1000:
            dataset = (name, X, y)
            datasets.append(dataset)

    return datasets

def get_datasets(typ):
    if 'Regression' in typ:
        return get_regression_datasets()
    elif 'Classification'in typ:
        return get_classification_datasets()

def get_regression_error(C, y):
    # Compute prediction error in the regression case
    def isclose(xs, ys):
        close = np.array([math.isclose(x, y) for x, y in zip(xs, ys)])
        return close

    close = np.logical_or(isclose(C[0], y), isclose(C[1], y))
    contain = np.logical_and(C[0] <= y, y <= C[1])

    error = 1. - np.logical_or(close, contain)
    return error

def get_classification_error(C, y):
    # Compute prediction error in the classification case
    n = C.shape[0]
    error = (1. - C[np.arange(n), y])
    return error

def get_error(typ, C, y):
    if 'Regression' in typ:
        return get_regression_error(C, y)
    elif 'Classification'in typ:
        return get_classification_error(C, y)

def get_regression_size(C):
    # Compute prediction set size in the regression case
    size = (C[1] - C[0])
    return size

def get_classification_size(C):
    # Compute prediction set size in the classification case
    size = C.sum(1)
    return size

def get_size(typ, C):
    if 'Regression' in typ:
        return get_regression_size(C)
    elif 'Classification'in typ:
        return get_classification_size(C)

def get_size_error(size, size_lower, size_upper):
    # Compute prediction set size interval error
    def isclose(xs, y):
        close = np.array([math.isclose(x, y) for x in xs])
        return close

    close = np.logical_or(isclose(size_lower, size), isclose(size_upper, size))
    contain = np.logical_and(size_lower <= size, size <= size_upper)

    size_error = 1. - np.logical_or(close, contain)
    return size_error

def get_mean_std(x):
    # Compute the mean and the standard deviation
    mean = x.mean()
    std = x.std(ddof=1)
    return mean, std


utils.py
Guneet-Dhillon/expected-conformal-prediction-set-size