Create Prediction Intervals such as Quantile Regression, Conformalized Quantile Regression, and bootstrapping using XGBoost
Project description
Based off the work done here:
- CQR: https://arxiv.org/abs/1905.03222
- CQR: poster concise explanation: https://github.com/yromano/cqr/blob/master/poster/CQR_Poster.pdf
macOS users install libomp before running the project:
brew install libomp
```python
%pip install xgboost_prediction_interval
import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
from xgboost_prediction_interval import XGBoostCQR, XGBoostQuantileRegressor, XGBoostBootstrap
Demo of how to use package:
Creating dummy data
np.random.seed(42)
n_samples = 1000
# Generate some random features
X1 = np.random.normal(loc=0, scale=1, size=n_samples)
X2 = np.random.normal(loc=5, scale=2, size=n_samples)
X3 = np.random.normal(loc=5, scale=2, size=n_samples)
X4 = np.random.normal(loc=5, scale=2, size=n_samples)
X5 = np.random.normal(loc=5, scale=2, size=n_samples)
X6 = np.random.normal(loc=5, scale=2, size=n_samples)
X7 = np.random.normal(loc=5, scale=2, size=n_samples)
# Create a target variable with some noise
# assume a linear relationship with some noise
y = 3 * X1 + 2 * X2 + np.random.normal(loc=0, scale=1, size=n_samples) +X3+X4+X5+X6+X7
df = pd.DataFrame({
"Feature1": X1,
"Feature2": X2,
"Feature3": X3,
"Feature4": X4,
"Feature5": X5,
"Feature6": X6,
"Feature7": X7,
"Target": y
})
print(df.head())
print(df.shape)
Feature1 Feature2 Feature3 Feature4 Feature5 Feature6 Feature7 \
0 0.496714 7.798711 3.649643 1.184385 3.273013 4.152481 2.771837
1 -0.138264 6.849267 4.710963 3.279230 4.937593 4.093172 3.738138
2 0.647689 5.119261 3.415160 4.172789 5.036034 1.408714 3.115880
3 1.523030 3.706126 4.384077 8.775375 5.945261 4.339820 3.904008
4 -0.234153 6.396447 1.212771 6.113106 2.266283 6.465658 4.571699
Target
0 32.904108
1 32.265157
2 30.044909
3 39.096159
4 33.427409
(1000, 8)
To import data saved locally uncomment the following and add full path to dataset in place of path_to_file...
# # Importing data saved locally:
# path = r"path_to_file..."
# data = pd.read_csv(path)
Splitting data
X_data = df.copy()
y_data = X_data.pop("Target")
X_train, X_test, y_train, y_test = train_test_split(X_data, y_data, test_size=0.25)
X_calibration, X_test, y_calibration, y_test = train_test_split(X_test, y_test, test_size=0.5)
xgboost_params = {
"max_depth": 4,
"eta": 0.05,
"gamma": 0.1,
"subsample": 0.8,
"min_child_weight": 100,
"colsample_bytree": 1,
"alpha": 5,
"lambda": 0,
"base_score": 0.57322098
}
cqr_xgboost = XGBoostCQR(
model_params=xgboost_params,
num_boost_round=int(100),
early_stopping_rounds=10,
alpha=0.9 # coverage probability for the conformal prediction part of CQR - should be set to your final desired coverage level
)
Using grid search to test different coverage probabilities for the quantile regression portion of CQR
Grid search is NOT a required step - can skip straight to fitting final PI model
The conformal prediction alpha ensures the PIs are padded to maintain the desired coverage probability even with a lower coverage specified for the quantile regression portion
lower_qr_quantile, upper_qr_quantile, conformity_score = cqr_xgboost.cqr_grid_search_alpha(
[0.1, 0.2, 0.4,], # lower quantile coverage probability levels
[0.6, 0.8, 0.9], # upper quantile coverage probability levels
X_train.values,
y_train.values,
X_calibration.values, # unseen holdout calibration set used for the conformal prediction portion
y_calibration.values,
X_test.values,
y_test.values,
0.25 # split size for the training data of the initial quantile regression
)
Evaluating --- Lower alpha: 0.1 --- Upper alpha 0.6
Evaluating --- Lower alpha: 0.1 --- Upper alpha 0.8
Evaluating --- Lower alpha: 0.1 --- Upper alpha 0.9
Evaluating --- Lower alpha: 0.2 --- Upper alpha 0.6
Evaluating --- Lower alpha: 0.2 --- Upper alpha 0.8
Evaluating --- Lower alpha: 0.2 --- Upper alpha 0.9
Evaluating --- Lower alpha: 0.4 --- Upper alpha 0.6
Evaluating --- Lower alpha: 0.4 --- Upper alpha 0.8
Evaluating --- Lower alpha: 0.4 --- Upper alpha 0.9
best CWC: 2.2498335501828367e-12
best coverage: 0.952
QR Alphas: [0.4, 0.6]
l_alpha u_alpha cwc coverage average_PI_width
0 0.1 0.6 3.123635e-12 0.960 19.704405
1 0.1 0.8 4.587595e-12 0.952 16.637457
2 0.1 0.9 7.390271e-12 0.936 16.710157
3 0.2 0.6 4.482768e-12 0.952 17.158796
4 0.2 0.8 3.692294e-12 0.960 16.109138
5 0.2 0.9 7.937115e-12 0.936 15.027287
6 0.4 0.6 2.249834e-12 0.976 16.465689
7 0.4 0.8 3.815250e-12 0.960 15.331766
8 0.4 0.9 4.629115e-12 0.952 16.430967
The output identifies the best quantile values for the quantile regression portion of the prediction interval (0.4 and 0.6)
- despite the quantile levels giving a prediction interval width of 0.6-0.4 = 0.2 the overall PI will still have a width of 0.9 due to the alpha defined at the start
Using the results of the grid search we can fit the CQR model with the optimum quantile levels
cqr_xgboost.fit(
X_train.values,
y_train.values,
X_calibration.values,
y_calibration.values,
0.25, # training/validation data split proportion used in the quantile regression portion
0.4, # lower quantile value for quantile regression portion
0.6, # lower quantile value for quantile regression portion
)
Evaluating --- Lower QR quantile: 0.4 --- Upper QR quantile 0.6 --- CQR alpha 0.9
cqr_preds = cqr_xgboost.predict(X_test)
cqr_preds
{'model_40': <xgboost.core.Booster object at 0x000001A4BE092630>, 'model_60': <xgboost.core.Booster object at 0x000001A4BE091F40>}
(np.float64(6.691801754053905), np.float64(6.691801754053905))
{'model_40_predictions': array([23.486837, 18.67393 , 32.484924, 27.112062, 28.803766, 27.276163,
22.676764, 31.016535, 33.065357, 29.88289 , 29.197748, 23.760366,
25.238892, 28.017939, 25.820185, 30.838179, 26.705332, 25.86778 ,
29.041285, 26.21456 , 30.407969, 27.181593, 29.95053 , 23.81866 ,
30.059622, 29.288862, 28.555567, 28.5591 , 28.319204, 21.84705 ,
24.585625, 26.863321, 25.627268, 24.067434, 26.742422, 26.236872,
33.069336, 24.5562 , 30.418066, 27.926493, 21.664396, 23.109285,
28.026175, 31.024427, 21.150793, 22.963375, 28.087214, 27.088976,
24.280598, 22.03742 , 22.731644, 31.532808, 22.613436, 20.83436 ,
29.681917, 29.690145, 28.256884, 30.230074, 31.142015, 29.283438,
26.716932, 25.042675, 23.072714, 19.717453, 24.86407 , 29.446177,
31.87668 , 26.152586, 29.279047, 29.858416, 26.830156, 25.162756,
27.63354 , 19.20571 , 24.69247 , 25.690577, 30.029097, 30.762701,
18.25805 , 27.989725, 24.016665, 26.895353, 30.669115, 29.7009 ,
30.946955, 32.54754 , 24.912683, 27.601076, 31.12282 , 29.557177,
24.888874, 25.8893 , 30.666895, 28.617975, 21.349886, 23.419872,
28.527033, 24.727304, 30.481913, 27.833437, 22.214617, 30.06967 ,
32.49518 , 26.203722, 23.8881 , 25.841269, 24.499422, 29.215921,
29.807978, 22.654224, 26.073137, 28.352068, 27.04648 , 30.973734,
31.945948, 26.06415 , 27.427275, 24.810951, 24.995848, 28.251188,
30.708555, 18.834095, 25.407259, 23.607534, 25.937593],
dtype=float32),
'model_60_predictions': array([38.722267, 33.717125, 46.950207, 42.872467, 44.487522, 43.517685,
36.2037 , 47.33104 , 49.405422, 45.57175 , 45.849594, 40.34343 ,
42.515392, 45.768757, 41.283974, 47.002872, 43.49001 , 40.70927 ,
44.985947, 42.033783, 45.824375, 43.750404, 45.824776, 38.606483,
44.9471 , 44.982548, 43.323177, 44.30441 , 46.64857 , 36.690807,
42.08947 , 43.14702 , 40.105347, 42.03816 , 42.343517, 42.82062 ,
47.231827, 39.99295 , 46.01658 , 44.035904, 40.032692, 39.847237,
43.620773, 45.87897 , 38.059414, 39.671467, 44.383564, 43.294315,
42.21498 , 37.707195, 39.172314, 48.78219 , 40.33695 , 35.455032,
45.304146, 44.3602 , 44.47025 , 44.866173, 47.033566, 45.58766 ,
43.899357, 39.882862, 38.573433, 36.740124, 42.331272, 44.73848 ,
48.787716, 42.771965, 43.78315 , 46.3867 , 43.134953, 38.920753,
42.35099 , 34.952 , 40.73547 , 41.644997, 45.818905, 48.361866,
34.787315, 42.48722 , 39.872192, 43.96925 , 46.86797 , 45.07383 ,
48.388138, 47.090015, 40.41863 , 44.758327, 45.3852 , 44.541466,
42.1466 , 43.65617 , 47.54282 , 44.77244 , 36.969772, 41.54296 ,
43.831978, 39.787712, 46.23097 , 44.46332 , 39.139694, 44.293724,
49.27917 , 42.516434, 39.81532 , 42.348747, 41.600716, 44.78814 ,
43.9265 , 38.41995 , 42.405285, 42.678352, 40.930042, 47.90765 ,
45.33097 , 41.709576, 40.599304, 40.70203 , 40.479744, 41.99648 ,
48.51149 , 36.695763, 41.383724, 39.570942, 41.64385 ],
dtype=float32)}
creates a dictionary for the upper and lower prediction interval with the following naming convention: model_{quantile_level}_predictions
- the key value changes respect to the quantile level specified during fit()
pd.DataFrame(cqr_preds).head()
.dataframe tbody tr th {
vertical-align: top;
}
.dataframe thead th {
text-align: right;
}
</style>
| model_40_predictions | model_60_predictions | |
|---|---|---|
| 0 | 23.486837 | 38.722267 |
| 1 | 18.673929 | 33.717125 |
| 2 | 32.484924 | 46.950207 |
| 3 | 27.112062 | 42.872467 |
| 4 | 28.803766 | 44.487522 |
lower_quantile_preds = cqr_preds["model_40_predictions"]
upper_quantile_preds = cqr_preds["model_60_predictions"]
cqr_xgboost.plot_pi_line_graph(y_test, lower_quantile_preds, upper_quantile_preds)
print(f"symmetric coverage {cqr_xgboost.coverage(y_test, lower_quantile_preds, upper_quantile_preds)}")
print(f"symmetric average_width {cqr_xgboost.average_width(lower_quantile_preds, upper_quantile_preds)}")
print(f"symmetric cwc {cqr_xgboost.cwc(0.9, y_test, lower_quantile_preds, upper_quantile_preds)}")
symmetric coverage 0.96
symmetric average_width 15.935728073120117
symmetric cwc 3.719721652657278e-12
Also have the option for an asymmetric loss function which optimises the upper and lower quantiles separately
# instantiate the model
asymmetric_cqr_xgboost = XGBoostCQR(model_params=xgboost_params,num_boost_round=int(100), alpha=0.90)
# grid search upper and lower quantiles
lower_qr_quantile, upper_qr_quantile, conformity_score = cqr_xgboost.cqr_grid_search_alpha(
[0.1, 0.2, 0.4,], # lower quantile coverage probability levels
[0.6, 0.8, 0.9], # upper quantile coverage probability levels
X_train.values,
y_train.values,
X_calibration.values, # unseen holdout calibration set used for the conformal prediction portion
y_calibration.values,
X_test.values,
y_test.values,
0.25 # split size for the training data of the initial quantile regression
)
Evaluating --- Lower alpha: 0.1 --- Upper alpha 0.6
Evaluating --- Lower alpha: 0.1 --- Upper alpha 0.8
Evaluating --- Lower alpha: 0.1 --- Upper alpha 0.9
Evaluating --- Lower alpha: 0.2 --- Upper alpha 0.6
Evaluating --- Lower alpha: 0.2 --- Upper alpha 0.8
Evaluating --- Lower alpha: 0.2 --- Upper alpha 0.9
Evaluating --- Lower alpha: 0.4 --- Upper alpha 0.6
Evaluating --- Lower alpha: 0.4 --- Upper alpha 0.8
Evaluating --- Lower alpha: 0.4 --- Upper alpha 0.9
best CWC: 2.2979702687205084e-12
best coverage: 0.944
QR Alphas: [0.1, 0.6]
l_alpha u_alpha cwc coverage average_PI_width
0 0.1 0.6 2.297970e-12 0.968 20.983686
1 0.1 0.8 5.785209e-12 0.944 16.820417
2 0.1 0.9 7.639414e-12 0.936 15.943439
3 0.2 0.6 5.581861e-12 0.944 17.615952
4 0.2 0.8 3.708949e-12 0.960 16.003836
5 0.2 0.9 6.238370e-12 0.944 15.047573
6 0.4 0.6 4.770299e-12 0.952 15.728806
7 0.4 0.8 3.616456e-12 0.960 16.588610
8 0.4 0.9 6.010613e-12 0.944 15.938600
The Output for Grid Search shows key metrics (Coverage-Width Criterion, Coverage, average PI width) - best alphas recommended based on lowest CWC which changes based on arbitrary choice of eta parameter in CWC formula (higher eta favours high coverage => wider PIs, lower eta favours narrow PIs => lower coverage) thus for best results experiment with different values of eta or use the full output to also compare coverage and PI width when settling on lower and upper quantile levels. The balance of coverage and PI width will vary based on application
# fit the model
asymmetric_cqr_xgboost.fit(
X_train.values,
y_train.values,
X_calibration.values,
y_calibration.values, 0.25,
lower_qr_quantile=0.4,
upper_qr_quantile=0.6,
conformity_score_method="asymmetric"
)
# make predictions
asym_cqr_preds = asymmetric_cqr_xgboost.predict(X_test)
# access the lower and upper prediction interval values for the X_test data
asym_lower_quantile_preds = asym_cqr_preds["model_40_predictions"]
asym_upper_quantile_preds = asym_cqr_preds["model_60_predictions"]
# plot the prediction interval against the actual y_test values
asymmetric_cqr_xgboost.plot_pi_line_graph(y_test, asym_lower_quantile_preds, asym_upper_quantile_preds)
# print out some key metrics evaluating the prediction interval
print(f"asymmetric coverage {asymmetric_cqr_xgboost.coverage(y_test, asym_lower_quantile_preds, asym_upper_quantile_preds)}")
print(f"asymmetric average_width {asymmetric_cqr_xgboost.average_width(asym_lower_quantile_preds, asym_upper_quantile_preds)}")
print(f"asymmetric cwc {asymmetric_cqr_xgboost.cwc(0.9, y_test, asym_lower_quantile_preds, asym_upper_quantile_preds)}")
Evaluating --- Lower QR quantile: 0.4 --- Upper QR quantile 0.6 --- CQR alpha 0.9
{'model_40': <xgboost.core.Booster object at 0x000001A4DA37B290>, 'model_60': <xgboost.core.Booster object at 0x000001A4DA3794F0>}
(np.float64(6.199072327530171), np.float64(7.433534967125333))
asymmetric coverage 0.96
asymmetric average_width 15.902301788330078
asymmetric cwc 3.725008720146514e-12
Results in a marginally lower coverage probability but a narrower average PI width
XGBoost Quantile Regression Demo
qr_xgboost = XGBoostQuantileRegressor(
model_params=xgboost_params,
num_boost_round=int(100),
early_stopping_rounds=10,
quantiles=[0.05, 0.95] # coverage probability for the conformal prediction part of CQR - should be set to your final desired coverage level
)
qr_xgboost.fit(X_train, y_train, validation_size=0.25)
qr_preds = qr_xgboost.predict(X_test)
# access the lower and upper prediction interval values for the X_test data
lower_qr_preds = qr_preds["model_5_predictions"]
upper_qr_preds = qr_preds["model_95_predictions"]
# plot the prediction interval against the actual y_test values
asymmetric_cqr_xgboost.plot_pi_line_graph(y_test, lower_qr_preds, upper_qr_preds)
# print out some key metrics evaluating the prediction interval
print(f"asymmetric coverage {asymmetric_cqr_xgboost.coverage(y_test, lower_qr_preds, upper_qr_preds)}")
print(f"asymmetric average_width {asymmetric_cqr_xgboost.average_width(lower_qr_preds, upper_qr_preds)}")
print(f"asymmetric cwc {asymmetric_cqr_xgboost.cwc(0.9, y_test, lower_qr_preds, upper_qr_preds)}")
asymmetric coverage 0.888
asymmetric average_width 14.869775772094727
asymmetric cwc 3.3716176135545774e-11
bootstrap = XGBoostBootstrap(
model_params={},
num_boost_round=int(100),
alpha=0.98
)
bootstrap.fit(X_train, y_train, n_bootstrap=100, sample_size_ratio=0.7)
bs_preds = bootstrap.predict(X_test)
# access the lower and upper prediction interval values for the X_test data
lower_qr_preds = bs_preds["model_1_predictions"]
upper_qr_preds = bs_preds["model_99_predictions"]
# plot the prediction interval against the actual y_test values
bootstrap.plot_pi_line_graph(y_test, lower_qr_preds, upper_qr_preds)
# print out some key metrics evaluating the prediction interval
print(f"Bootstrap coverage {bootstrap.coverage(y_test, lower_qr_preds, upper_qr_preds)}")
print(f"Bootstrap average_width {bootstrap.average_width(lower_qr_preds, upper_qr_preds)}")
print(f"Bootstrap cwc {bootstrap.cwc(0.9, y_test, lower_qr_preds, upper_qr_preds)}")
----- Training model 1 / 100 -----
----- Training model 2 / 100 -----
----- Training model 3 / 100 -----
----- Training model 4 / 100 -----
----- Training model 5 / 100 -----
----- Training model 6 / 100 -----
----- Training model 7 / 100 -----
----- Training model 8 / 100 -----
----- Training model 9 / 100 -----
----- Training model 10 / 100 -----
----- Training model 11 / 100 -----
----- Training model 12 / 100 -----
----- Training model 13 / 100 -----
----- Training model 14 / 100 -----
----- Training model 15 / 100 -----
----- Training model 16 / 100 -----
----- Training model 17 / 100 -----
----- Training model 18 / 100 -----
----- Training model 19 / 100 -----
----- Training model 20 / 100 -----
----- Training model 21 / 100 -----
----- Training model 22 / 100 -----
----- Training model 23 / 100 -----
----- Training model 24 / 100 -----
----- Training model 25 / 100 -----
----- Training model 26 / 100 -----
----- Training model 27 / 100 -----
----- Training model 28 / 100 -----
----- Training model 29 / 100 -----
----- Training model 30 / 100 -----
----- Training model 31 / 100 -----
----- Training model 32 / 100 -----
----- Training model 33 / 100 -----
----- Training model 34 / 100 -----
----- Training model 35 / 100 -----
----- Training model 36 / 100 -----
----- Training model 37 / 100 -----
----- Training model 38 / 100 -----
----- Training model 39 / 100 -----
----- Training model 40 / 100 -----
----- Training model 41 / 100 -----
----- Training model 42 / 100 -----
----- Training model 43 / 100 -----
----- Training model 44 / 100 -----
----- Training model 45 / 100 -----
----- Training model 46 / 100 -----
----- Training model 47 / 100 -----
----- Training model 48 / 100 -----
----- Training model 49 / 100 -----
----- Training model 50 / 100 -----
----- Training model 51 / 100 -----
----- Training model 52 / 100 -----
----- Training model 53 / 100 -----
----- Training model 54 / 100 -----
----- Training model 55 / 100 -----
----- Training model 56 / 100 -----
----- Training model 57 / 100 -----
----- Training model 58 / 100 -----
----- Training model 59 / 100 -----
----- Training model 60 / 100 -----
----- Training model 61 / 100 -----
----- Training model 62 / 100 -----
----- Training model 63 / 100 -----
----- Training model 64 / 100 -----
----- Training model 65 / 100 -----
----- Training model 66 / 100 -----
----- Training model 67 / 100 -----
----- Training model 68 / 100 -----
----- Training model 69 / 100 -----
----- Training model 70 / 100 -----
----- Training model 71 / 100 -----
----- Training model 72 / 100 -----
----- Training model 73 / 100 -----
----- Training model 74 / 100 -----
----- Training model 75 / 100 -----
----- Training model 76 / 100 -----
----- Training model 77 / 100 -----
----- Training model 78 / 100 -----
----- Training model 79 / 100 -----
----- Training model 80 / 100 -----
----- Training model 81 / 100 -----
----- Training model 82 / 100 -----
----- Training model 83 / 100 -----
----- Training model 84 / 100 -----
----- Training model 85 / 100 -----
----- Training model 86 / 100 -----
----- Training model 87 / 100 -----
----- Training model 88 / 100 -----
----- Training model 89 / 100 -----
----- Training model 90 / 100 -----
----- Training model 91 / 100 -----
----- Training model 92 / 100 -----
----- Training model 93 / 100 -----
----- Training model 94 / 100 -----
----- Training model 95 / 100 -----
----- Training model 96 / 100 -----
----- Training model 97 / 100 -----
----- Training model 98 / 100 -----
----- Training model 99 / 100 -----
----- Training model 100 / 100 -----
Bootstrap coverage 0.896
Bootstrap average_width 7.640153408050537
Bootstrap cwc 3.4321849985369016e-11
Bootstrap undercovers with the default parameters - introducing randomness via XGBoosts hyperparameters (subsample, colsample_bytree, colsample_bylevel , ...) may reduce this effect
bootstrap.plot_pi_line_graph(y_test, lower_qr_preds, upper_qr_preds, X_test["Feature1"])
bootstrap.plot_coverage_probability_binned(y_test, lower_qr_preds, upper_qr_preds, x=X_test["Feature1"])
bootstrap.plot_pi_width(lower_qr_preds, upper_qr_preds, x=X_test["Feature1"])
Project details
Release history Release notifications | RSS feed
Download files
Download the file for your platform. If you're not sure which to choose, learn more about installing packages.
Source Distributions
Built Distribution
Filter files by name, interpreter, ABI, and platform.
If you're not sure about the file name format, learn more about wheel file names.
Copy a direct link to the current filters
File details
Details for the file prediction_interval-0.1.0-py3-none-any.whl.
File metadata
- Download URL: prediction_interval-0.1.0-py3-none-any.whl
- Upload date:
- Size: 9.5 kB
- Tags: Python 3
- Uploaded using Trusted Publishing? No
- Uploaded via: twine/6.1.0 CPython/3.13.2
File hashes
| Algorithm | Hash digest | |
|---|---|---|
| SHA256 |
6474a985bf725531f2e27d88f60b0d3a4db1b99e1c612c094098b361ba726fc5
|
|
| MD5 |
d87061dc0829e0e6d995d54efeef275e
|
|
| BLAKE2b-256 |
756553be116f1899215f9116f54a906b2356de1e883172b3015782fa9d62af63
|