  Let $\mathcal{P}_{\beta}^{(V)} (N_{\cal I})$ be the probability that a
$N\times N$ $\beta$-ensemble of random matrices with confining potential $V(x)$
has $N_{\cal I}$ eigenvalues inside an interval ${\cal I}=[a,b]$ of the real
line. We introduce a general formalism, based on the Coulomb gas technique and
the resolvent method, to compute analytically $\mathcal{P}_{\beta}^{(V)}
(N_{\cal I})$ for large $N$. We show that this probability scales for large $N$
as $\mathcal{P}_{\beta}^{(V)} (N_{\cal I})\approx \exp\left(-\beta N^2
\psi^{(V)}(N_{\cal I} /N)\right)$, where $\beta$ is the Dyson index of the
ensemble. The rate function $\psi^{(V)}(k_{\cal I})$, independent of $\beta$,
is computed in terms of single integrals that can be easily evaluated
numerically. The general formalism is then applied to the classical
$\beta$-Gaussian (${\cal I}=[-L,L]$), $\beta$-Wishart (${\cal I}=[1,L]$) and
$\beta$-Cauchy (${\cal I}=[-L,L]$) ensembles. Expanding the rate function
around its minimum, we find that generically the number variance ${\rm
Var}(N_{\cal I})$ exhibits a non-monotonic behavior as a function of the size
of the interval, with a maximum that can be precisely characterized. These
analytical results, corroborated by numerical simulations, provide the full
counting statistics of many systems where random matrix models apply. In
particular, we present results for the full counting statistics of zero
temperature one-dimensional spinless fermions in a harmonic trap.
