XGBoost Classifier not capturing extreme probabilities I'm using XGBoost for a binary classification task—trying to predict whether team A will beat team B given the score of the game and the time left. I know for certain score-time combinations, the probability of a success (team A wins) is essentially 0 or 1 (e.g. up 20 points with a minute left, it should be ~0.999). This is true empirically (in previous games) regardless of the other features (team quality, pace of play, etc.)
The problem: When I look at fitted values given these score-time combos, the output tends to be around 0.03 or 0.97, rather than .0001 or .9999. I've tried loads of different combinations of parameters and can't get the model to output something close to 0 or 1.
A few more details:

*

*I have about 630k observations, about 55% of which are successes. I also have at least tens of thousands where I believe the probability should be > .99 or < .01.

*Parameters I've messed with + the ranges I've tried: max_depth (from 4 to 12), n_estimators (up to a couple hundred), eta (.001 to .3), min_child_weight (up to a couple hundred), reg_lambda (1-5), gamma (0-5), colsample_byX (0.7-1), and subsample (0.7-1)

*I have not messed with scale_pos_weight since my understanding is that this parameter is for imbalanced datasets. My dataset isn't imbalanced overall, it's just imbalanced at certain places in the feature space

*Loss function is binary:logistic
Any thoughts? Thanks a ton!

As @Sycorax and @BenReiniger pointed out, the problem is that the probabilities are not calibrated (or not as calibrated as well as you'd prefer). Here is how you could calibrate the XGBoost probabilities. Use the following model:
P(y|x) = 1/(1+exp(-(a+x)))
where x is the logit function of the original probabilities produced by XGBoost:
logit = log(p/(1-p))
and y are the same outcomes you are already using. This is based on the paper by van den Goorbergh et al., "The harm of class imbalance corrections for risk prediction models: illustration and simulation using logistic regression", arXiv 2022 (see the Methods section). In my experience works well.
You can implement this in R using this statement:
recal_mod = glm(y ~ 1, offset = logit, family = "binomial")
Note that this model is a logistic regression whereas the independent variable has a fixed weight of 1 and only the intercept is fit. To my knowledge this type of model is not supported in python scikit-learn
LogisticRegression() function, so I use R. Note also that it uses lot of RAM for large datasets, so you may want to to downsample your data.
The following graphs illustrate how it works. The first graph is a calibration curve before recalibration. The second graph is the calibration curve after recalibration using the above model. I don't know if it will achieve the calibration that you want but it's worth trying


