Finding largest subset that matches moments
I tried formulating this as an integer linear program, but ran into problems, because the variance is not a linear function: it is quadratic.



Here was my attempt.

Introduce variables $x_b$ (one for each $b \in B$), such that $x_b=1$ if $b \in C$ and $x_b=0$ if $b \notin C$.  Thus, we have $|B|$ new zero-or-one variables.  We can force them to be zero or one by adding the constraints $0 \le x_b \le 1$ for each $b \in B$.

Next, we'll force the mean of $C$ to be in the accepted range.  Introduce a new unknown $\mu$ and add the constraint

$$\mu = (\sum_{b \in B} x_b b)/|B|.$$

This forces $\mu$ to be the mean of $C$.  Also, it is expressible in a linear equation, since everything except for $\mu$ and the $x$'s are constants.  Finally, add the two linear inequalities

$$\text{mean}(A)-\delta \le \mu \le \text{mean}(A)+\delta.$$

But it's not clear how to force the standard deviation of $C$ to be in the accepted range.  Let $\sigma = \text{std}(C)$, i.e., $\sigma$ denotes the standard deviation of $C$.  We could try to use the fact that

$$\sigma^2 =  \frac{1}{N} \sum_{c \in C} c^2 - \frac{1}{N^2} \left( \sum_{c \in C} c \right)^2,$$

but now when we substitute in with our definitions of $x_b$, we get something ugly:

$$\sigma^2 =  \frac{1}{N} \sum_{b \in B} b^2 x_b - \frac{1}{N^2} \left( \sum_{b \in B} b x_b \right)^2,$$

where the right-hand term has something quadratic.  Worse still, $N$ also depends upon the $x$'s, which takes this even further away from something linear:

$$S = \sum_{b \in B} x_b.$$

Alternatively, we could try to use the fact that

$$\sigma^2 = \frac{1}{N} \sum_{c \in C} (c - \mu)^2,$$

but when we substitute in, again we get something that's not linear:

$$\sigma^2 = \frac{1}{N} \sum_{b \in B} x_b (b - \mu)^2,$$

where the non-linearity comes from the fact that $N$ is itself a function of the $b$'s, as well as the fact that the right-hand side involves a product of two linear variables ($x_b$ and $\mu$).

We want

$$\text{std}(A) - \epsilon \le \sigma^2 \le \text{std}(A) + \epsilon.$$

Plugging in, we find that we need

$$(\text{std}(A) - \epsilon)N \le \sum_{b \in B} x_b (b - \mu)^2 \le (\text{std}(A) + \epsilon)N.$$

This is not linear, because $x_b (b - \mu)^2$ is not a linear function of $x_b,\mu$.

We could deal with the fact that $N$ is not a constant by enumerating all candidate values of $N$ and solving an ILP for each.  However, that still doesn't deal with the problem that $x_b (b - \mu)^2$ is non-linear.

So, I'm stuck.  I wasn't able to formulate this cleanly as an ILP.