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.