Reproducing results of the paper: Distributional regression#
This notebook contains the code used to reproduce the results presented in our paper comparing different methods for distributional regression.
[2]:
import numpy as np
np.random.seed(1)
import os
import sys
import matplotlib.pyplot as plt
# Get the directory where the current notebook is located
NOTEBOOK_DIR = os.getcwd()
# Go up two levels (adjust the '..' count as needed)
REPO_ROOT = os.path.abspath(os.path.join(NOTEBOOK_DIR, '../../../../'))
# Add to sys.path if not already there
if REPO_ROOT not in sys.path:
sys.path.insert(0, REPO_ROOT)
from DisTreebution.UQ.UQ import UQ
import pickle
import random
import argparse
from DisTreebution.UQ.utils import compute_crps,load_dataset
1. Computing the CRPS metric for the following different methods: CRPS-RF, PMQRF and SKLearn#
[ ]:
def process_dataset(data_path, name_dataset):
for ite in range(100):
ntrain = 1000
if name_dataset=='red_wine':
ntest = 1000
ntest = 500
np.random.seed(ite)
random.seed(ite)
X, y = load_dataset(data_path, name_dataset)
N = X.shape[0]
N = min(N,ntrain+3000)
idxs = np.arange(len(y))
np.random.shuffle(idxs)
print(f'Processing {name_dataset}')
x_train = X[idxs[:N][:ntrain], :]
y_train = y[idxs[:N][:ntrain]]
x_test = X[idxs[:N][(ntrain):], :]
y_test = y[idxs[:N][(ntrain):]]
nTrees = 100
alpha = 0.1
ls_q = [0.05 * i for i in range(1, 20)]
results = []
# Baseline SKLearn Quantile Regressor
from sklearn.linear_model import QuantileRegressor
marginal_level = np.zeros(len(ls_q))
sample2quantiles = np.zeros((len(y_test), len(ls_q)))
for i_q, q in enumerate(ls_q):
model = QuantileRegressor(quantile=q, alpha=0)
model.fit(x_train, y_train)
predictions = model.predict(x_test)
sample2quantiles[:,i_q] = predictions
crps_value = compute_crps(sample2quantiles, y_test, ls_q)
results.append({
'dataset': name_dataset,
'type_tree': 'SKLearn',
'type_conformal': None,
'type_aggregation_trees': None,
'nested_set': None,
'groupconf': None,
'IG_biais_correction': None,
'split': None,
'crps_value': crps_value,
'ntrain': ntrain,
'alpha':alpha
})
model_configs = []
for type_tree in ['CRPS','PMQRT']:
for type_aggregation_trees in ['vr']:
for split in [False, True]:
model_configs.append(
{
'type_tree': type_tree,
'type_conformal': None,
'type_aggregation_trees': type_aggregation_trees,
'nested_set': None,
'groupconf': False,
"IG_biais_correction": "LOO",
"split":split
}
)
nested_set = None
groupconf = False
for config in model_configs:
type_tree, type_conformal = config['type_tree'], config['type_conformal']
type_aggregation_trees = config['type_aggregation_trees']
IG_biais_correction = config["IG_biais_correction"]
split = config["split"]
if split and type_tree!="PMQRT":
pass
else:
print(type_tree, type_conformal, nested_set, "IG_biais_correction", IG_biais_correction)
params = {'nTrees': nTrees, 'max_depth': 10, 'min_samples_split': 10,
'IG_biais_correction':IG_biais_correction, 'list_distri_low_quantiles': [0.01 * i for i in range(1, 20)]}
treeID2quantiles_train = None
if type_tree == "PMQRT":
if not (split):
treeID2quantiles_train = {ID:[0.05 * i for i in range(1, 20)] for ID in range(nTrees)}
else:
treeID2quantiles_train = {ID: [0.05 * i for i in range(1, 10)] for ID in range(nTrees//2)}
treeID2quantiles_train.update({ID: [0.05 * i for i in range(10,20)] for ID in range(nTrees//2,nTrees)})
params.update({'treeID2quantiles_train': treeID2quantiles_train})
model = UQ(type_tree=type_tree, nested_set=nested_set, type_conformal=type_conformal, group_coverage=groupconf, type_aggregation_trees=type_aggregation_trees, params=params)
print('Training')
trees, sample2calib_trees = model.train_trees(x_train, y_train)
print('Inference')
sample2quantiles = model.get_quantile_estimate(trees, x_test, quantiles=ls_q)
crps_value = compute_crps(sample2quantiles, y_test, ls_q)
results.append({
'dataset': name_dataset,
'type_tree': type_tree,
'type_conformal': type_conformal,
'type_aggregation_trees': type_aggregation_trees,
'nested_set': nested_set,
'groupconf': groupconf,
'IG_biais_correction': IG_biais_correction,
'split': split,
'crps_value': crps_value,
'ntrain': ntrain,
'alpha':alpha
})
# Save all
import pandas as pd
results_df = pd.DataFrame(results)
#results_df.to_csv(os.path.join(save_path, f'results_final_crps_{name_dataset}_ite_{ite}.csv'), index=False)
print(f"Saved {len(results)} results to results_{name_dataset}_{ite}.csv")
if __name__ == '__main__':
parser = argparse.ArgumentParser(description="Configure regression tree and conformalization options.")
parser.add_argument("--datasetID", type=int, default=None,
help="datasetID")
parser.add_argument("--data_path", type=int, default=None,
help="data path")
args = parser.parse_args()
datasets = ["abalone", "gpu", "gas_turbine", "combined_cycle_power_plant", "red_wine", "white_wine"]
name_dataset = datasets[args.datasetID]
data_path = args.data_path
process_dataset(data_path, name_dataset)
2. Computing the CRPS metric for QRF#
[ ]:
def process_dataset(name_dataset):
for ite in range(100):
ntrain = 1000
if name_dataset=='red_wine':
ntest = 1000
ntest = 500
np.random.seed(ite)
random.seed(ite)
X, y = load_dataset(name_dataset)
N = X.shape[0]
N = min(N,ntrain+3000)
idxs = np.arange(len(y))
np.random.shuffle(idxs)
print(f'Processing {name_dataset}')
x_train = X[idxs[:N][:ntrain], :]
y_train = y[idxs[:N][:ntrain]]
x_test = X[idxs[:N][(ntrain):], :]
y_test = y[idxs[:N][(ntrain):]]
nTrees = 100
alpha = 0.1
ls_q = [0.05 * i for i in range(1, 20)]
results = []
from sklearn.linear_model import QuantileRegressor
marginal_level = np.zeros(len(ls_q))
sample2quantiles = np.zeros((len(y_test), len(ls_q)))
for i_q, q in enumerate(ls_q):
model = QuantileRegressor(quantile=q, alpha=0)
model.fit(x_train, y_train)
predictions = model.predict(x_test)
sample2quantiles[:,i_q] = predictions
crps_value = compute_crps(sample2quantiles, y_test, ls_q)
results.append({
'dataset': name_dataset,
'type_tree': 'SKLearn',
'type_conformal': None,
'type_aggregation_trees': None,
'nested_set': None,
'groupconf': None,
'IG_biais_correction': None,
'split': None,
'crps_value': crps_value,
'ntrain': ntrain,
'alpha':alpha
})
model_configs = []
for type_tree in ['RT']:
for type_aggregation_trees in ['vr','vr-avg']:
model_configs.append(
{
'type_tree': type_tree,
'type_conformal': None,
'type_aggregation_trees': type_aggregation_trees,
'nested_set': None,
'groupconf': False,
"IG_biais_correction": None,
"split":False
}
)
nested_set = None
groupconf = False
for config in model_configs:
type_tree, type_conformal = config['type_tree'], config['type_conformal']
type_aggregation_trees = config['type_aggregation_trees']
IG_biais_correction = config["IG_biais_correction"]
split = config["split"]
print(type_tree, type_conformal, nested_set, "IG_biais_correction", IG_biais_correction)
params = {'nTrees': nTrees, 'max_depth': 10, 'min_samples_split': 10,
'IG_biais_correction':IG_biais_correction, 'list_distri_low_quantiles': [0.01 * i for i in range(1, 20)]}
model = UQ(type_tree=type_tree, nested_set=nested_set, type_conformal=type_conformal, group_coverage=groupconf, type_aggregation_trees=type_aggregation_trees, params=params)
print('Training')
trees, sample2calib_trees = model.train_trees(x_train, y_train)
print('Inference')
sample2quantiles = model.get_quantile_estimate(trees, x_test, quantiles=ls_q)
crps_value = compute_crps(sample2quantiles, y_test, ls_q)
results.append({
'dataset': name_dataset,
'type_tree': type_tree,
'type_conformal': type_conformal,
'type_aggregation_trees': type_aggregation_trees,
'nested_set': nested_set,
'groupconf': groupconf,
'IG_biais_correction': IG_biais_correction,
'split': split,
'crps_value': crps_value,
'ntrain': ntrain,
'alpha':alpha
})
# Save all
import pandas as pd
results_df = pd.DataFrame(results)
results_df.to_csv(os.path.join(save_path, f'results_final_qrf_{name_dataset}_ite_{ite}.csv'), index=False)
print(f"Saved {len(results)} results to results_{name_dataset}_{ite}.csv")
if __name__ == '__main__':
parser = argparse.ArgumentParser(description="Configure regression tree and conformalization options.")
parser.add_argument("--datasetID", type=int, default=None,
help="datasetID")
parser.add_argument("--data_path", type=int, default=None,
help="data path")
args = parser.parse_args()
datasets = ["abalone", "gpu", "gas_turbine", "combined_cycle_power_plant", "red_wine", "white_wine"]
name_dataset = datasets[args.datasetID]
data_path = args.data_path
process_dataset(data_path, name_dataset)
[ ]: