"Automatically" detect biased subsets in probability distribution Background:
Suppose we have a model generating probabilities conditional on a state vector - for simplicity we can just assume the outcome is 0 or 1 (imagine for example simple logistic regression):
$P(X = 1|S = s) = f(s)$.
An obvious sanity check is to make sure the model is "globally" unbiased by which I mean the following.
Suppose we have a set of datapoints $\{(S_i, X_i)\}_{i=1}^N$.
Then we want
$\sum_{i=1}^N \mathbb{1}_{\{X_i = 1\}} \approx \sum_{i=1}^N f(S_i)$.
But even if this holds, there could be subsets of $T \subset S$ for which it does not, i.e. conditional on the state residing in $T$ we have a bias:
$\sum_{\{i: s_i \in T\}} \mathbb{1}_{\{X_i = 1\}} \not\approx \sum_{\{i: s_i \in T\}} f(S_i)$.
Question:
What are the kinds of algorithms that apply to superficially healthy multiclass data that can identify problematically biased subsets?
I have some ideas for how one could go about this, but I wonder if there are some standard algorithms to find such subsets.
Update:
From the answers I got it seems a more concrete version of the problem would help.
Assume $S \in \mathbb{R}^k$ and our $f$ was obtained by logistic regression, i.e. for $s = \{s_j\}_{j=1}^{k}$
$f(s) = \frac{1}{1 + e^{-\sum_{j = 1}^k \beta_j s_j}  }$
where the $\beta_j$ are obtained by maximum likelihood as usual.
And I want to find subsets as natural as possible where the model has a significant bias and hence might be problematic.
Question 1: What do I mean by "significant" bias?
Answer 1: Given a subset $T \subset S$, we can view it as us having a null hypothesis that we have a collection of Bernoulli variables $\{ \text{Bern}(f(s_i)) \}_{\{i: s_i \in T\}}$
and we can compute in which percentile the observed number of variables having value 1 (call it $\text{prc}(T))$ resides given this null hypothesis - the more extreme the percentile, the more "significant" the bias.
Question 2: What do I mean by "natural" subsets?
Answer 2: Anything that is describable in terms of the state space really. Imagine e.g. that inside a hypersphere or radius $R$ we have a large bias, $T = \{ S : \sum_i s_i^2 \leq R^2 \}$.
My first idea was basically to use a greedy approach applied to regression/classification trees.
This would yield rectangular areas where the bias is potentially large.
But I'm curious of other approaches, and in particular which are the most commonly used for this type of - in my opinion - very natural problems.

Your proposal is quite different to what we usually do in binary regression problems
Your question uses a kind of binary regression model but it appears that your model has only a single regression function (as opposed to a parametric family of such functions).  At present it is also unclear exactly what you will use this proposed metric for, so it is difficult to give effective advice on the usefulness (or lack thereof) of your proposal.  In any case, I will take your question on its own terms and show how you can simplify this problem and identify the sets of interest to you, but I will also make some comments along the way comparing your approach to what we usually do in binary regression.  (Also, there is no need to assume that this is a logistic regression, so I will drop that assumption.)

Consider the case of a binary regression model with an assumed regression function $f$ (which may or may not be the logistic regression), giving $\mathbb{P}(X_i=1 | S_i=s) = \mathbb{E}(X_i | S_i=s) = f(s)$.  This regression model has sampling distribution given by:
$$p(X_1=x_1,...,X_n=x_n | \mathbf{S}=\mathbf{s}) = \prod_{i=1}^n \text{Bern}(x_i |f(s_i)).$$
Using the observed data points $(x_1,s_1),...,(x_n,s_n)$ we can easily compute the deviations of the actual response from its expected value as:
$$r_i \equiv x_i - f(s_i).$$
These give us some crude "residuals" for the model, measuring the deviation of the binary response variable from its expectation under the model.  (Since there is only a single regression function in this problem there is no estimated regression function; consequently, there is no distinction in this case between the "error terms" in the model and the "residuals".  You could just as easily call these quantities the error terms if you prefer, since there is no distinction in this case.)   You can use these quantities to simplify the expressions in your question.  To do this, observe for any set $\mathcal{T} \subseteq \{ 1,...,n \}$ (including the full set of values) we have:$^\dagger$
$$\sum_{i \in \mathcal{T}} r_i \approx 0
\quad \quad \quad \iff \quad \quad \quad 
\sum_{i \in \mathcal{T}} x_i \approx \sum_{i \in \mathcal{T}} f(s_i).$$

Your proposed analysis: In your question, you say that you want to have $\sum_{i \in \mathcal{T}} x_i \approx \sum_{i \in \mathcal{T}} f(s_i)$ as a test to see if the model is unbiased.  Since this is equivalent to $\sum_{i \in \mathcal{T}} r_i \approx 0$, you are effectively looking at whether the sum of residuals is near zero.  This is quite a crude way to look for bias in the model --- it is not wrong, but it is only going to detect a narrow set of forms of bias that are systematic across enough of the explanatory values that it manifests in large positive or negative totals for the residuals.  Your approach will not generally detect forms of bias where the regression is positively biased in some parts of the function and negatively biased in others (such that it "cancels out" in the residual sum).
Here it is worth noting what we usually do in regression analysis when we want to know if the posited regression function is reasonable or not.  What we usually do is to look at the residual plot showing the relationship between the residuals and the fitted values.  (In the case of a regression with a single explanatory variable, on can also look at the residuals against that explanatory variable.)  This type of plot will usually allow us to diagnose systematic bias in the regression function.  Moreover, it can be used to see if the regression function is biased in just some parts of the range of the explanatory variables, etc.  In the case of binary regression it is a bit tricky,  owing to the binary response, but the residual plot is still useful for this purpose if examined with care.
A secondary aspect of your approach that is potentially problematic is that you have not really specified how you will tell if the residual sum of interest is "close enough" to zero that you think it is unbiased.  To do this you would need to derive the distributional properties of the residuals, so that you can see the distribution of the residual sum (or at least some approximation to this distribution) under the model (i.e., assuming the specified regression function is the correct form).
Finally, while you have proposed identifying subsets of values where the residual sum deviates substantially from zero, you have not really described how you would apply this analysis.  As you correctly recognise in your question, even if the model is correct (and therefore unbiased) you can certainly cherry-pick subsets of values where the residual sum is substantially different from zero.  (I will show you how to do this below.)  Consequently, any attempt to use cherry-picked subsets like this to try to identify local bias in the regression is likely to be a very bad procedure; so it is unclear what the point is of identifying these subsets of values.

Identifying subsets of values of interest: Notwithstanding the above critique, I will take your question on its own terms and show you how you can compute the subsets of interest.  It is simple to use a sorting algorithm to obtain the ordered residuals, and it is also simple to identify where these ordered residuals cross zero.  Consequently, it is simple to compute index values $i_1,...,i_n$ and the value $k$ giving the following order:
$$r_{i_1} \leqslant \cdots \leqslant r_{i_k} < 0 \leqslant r_{i_{k+1}} \leqslant \cdots \leqslant r_{i_n}.$$
If you want to identify a subset of data points with an extreme negative residual sum, you can take any subset $\mathcal{T} = \{ i_1,...,i_r \}$ for some $r \leqslant k$.  Similarly, if you want to identify a subset of data points with an extreme positive residual sum, you can take any subset $\mathcal{T} = \{ i_{r+1},...,i_n\}$ for some $r \geqslant k$.  In both cases, you start by adding the data point with the most extreme negative/positive residual, then the second-most-extreme, and so on.
It is simple to construct an algorithm to do this.  I will illustrate this below in R.  First, I generate some mock data from a logistic regression model with a fixed regression function.  Then I show how to compute subsets giving extreme negative/positive residual sums.
#Generate logistic regression data with a fixed regression function
set.seed(1)
n  <- 100
S  <- rnorm(n)
fS <- 1/(1+exp(-S))
X  <- rep(0, n)
for (i in 1:n) { X[i] <- 1*(runif(1) <= fS[i]) }

#Compute the residuals/error terms
R  <- X - fS

#Order the residuals and count negative residuals
ORDER <- order(R)
k     <- sum(R<0)

#Identify subsets with extreme negative/positive residual sum
if (k > 0) { SUB.NEG <- sort(ORDER[1:k]) }     else { SUB.NEG <- integer(0) }
if (k < n) { SUB.POS <- sort(ORDER[(k+1):n]) } else { SUB.POS <- integer(0) }

#Show negative residual sum
sum(R[SUB.NEG])
[1] -18.3847

#Show positive residual sum
sum(R[SUB.POS])
[1] 25.9876


$^\dagger$ Note that the use of the indicator function for $x_i$ in your sums is unecessary, since $x_i$ is already a binary variable.  Consequently, this expression does not use indicator functions.
