Given a collection $\{λ_1, \dots, λ_n\} $ of real numbers, there is a canonical probability distribution on the set of real symmetric or complex Hermitian matrices with eigenvalues $λ_1,\ldots,λ_n$. In this paper, we study various features of random matrices with this distribution. Our main results show that under mild conditions, when $n$ is large, linear functionals of the entries of such random matrices have approximately Gaussian joint distributions. The results take the form of upper bounds on distances between multivariate distributions, which allows us also to consider the case when the number of linear functionals grows with $n$. In the context of quantum mechanics, these results can be viewed as describing the joint probability distribution of the expectation values of a family of observables on a quantum system in a random mixed state. Other applications are given to spectral distributions of submatrices, the classical invariant ensembles, and to a probabilistic counterpart of the Schur--Horn theorem, relating eigenvalues and diagonal entries of Hermitian matrices.