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.
