Probability that at least one sample lies within a given st.d.? I think my question is similar to this one, but different in that I consider a set of realisations, not only one. Sorry if this question is really easy, I'm just not sure how to go rigorously about it.
Say I have $N$ realisations from a multivariate normal distribution $\mathcal{N}(\mu,\Sigma)$. Intuitively, I would expect that the larger $N$, and the more likely I would be to get at least one realisation lying within $\epsilon$ st.d. of the mean $\mu$. I feel that this mixes discrete and continuous probabilities, and I need help to understand how to come up with what I guess is an expectation formula in this case. 
Note that while it is clear what "within $\epsilon$ st.d." means in the one-dimensional case, the multidimensional counterpart corresponds to "within the range $[0,\epsilon]$ of the Mahalanobis distance". 

If I had to have a go at it, I would find the probability of none being within, as the inverse cumulated density at $\epsilon$ st.d. (the probability of NOT being within), raised to the $N$, or:
$$
\left( \frac{1-\mathrm{erf}(\epsilon)}{2} \right)^N
$$
So the probability of at least one, would be one minus that, is that right?
 A: One-dimensional case
The probability of observing $x$ within the interval $\mu \pm \epsilon\sigma$ in a single draw is
$$ \Pr(X \in [\mu-\epsilon\sigma, \mu+\epsilon\sigma]) = \Pr(|X-\mu| \le \epsilon\sigma) $$
The event "at least one within this interval in $N$ draws" is the complement of "none within this interval in $N$ draws", so we first need to calculate the probability of not being in the interval:
$$ 1 - \Pr(|X-\mu| \le \epsilon\sigma) = \Pr\left( \left|\frac{X - \mu}{\sigma} \right| > \epsilon \right) = p $$
and then consider this happening $N$ times consecutively (but independently), which simply corresponds to the probability $p^N$ (cf formula in the OP). Finally, the event of interest is the complement of this, so $1-p^N$ is the answer.
As pointed out in comment by Dilip Sarwate, we have in the one-dimensional case
$$ \Pr\left( \left|\frac{X - \mu}{\sigma} \right| > \epsilon \right) = \Phi(\epsilon)-\Phi(-\epsilon) = 2\Phi(\epsilon)-1 = \mathrm{erf}(\epsilon/\sqrt{2}) $$
where we used the fact that $\Phi(-x) = 1-\Phi(x)$ for symmetric distributions.
Multi-dimensional case
The approach is similar, but now $p$ has a different formulation
$$ p=\Pr\left( \left|\frac{X - \mu}{\sigma} \right| > \epsilon \right) \quad\to\quad 
p=\Pr\left( \mathcal{M}_{\mu,\Sigma}(X)\in [0,\epsilon] \right) $$
using the Mahalanobis distance $ \mathcal{M}_{\mu,\Sigma}(x) = \sqrt{ (x-\mu)^t \Sigma^{-1} (x-\mu) } $. Furthermore the "cumulative" density function doesn't make sense anymore (there is no obvious partial order in dimension n), but I think the previous result $p = \mathrm{erf}(\epsilon/\sqrt{2})$ still holds regardless of the number of dimensions, although I don't know how to prove it.
