Here's one way you can think of it - we will specialize to a finite dimensional Hilbert space $\mathcal H$ for simplicity. Abstractly, the state of a quantum system should be identified with a probability distribution for every possible measurement one could make.
A given self-adjoint operator $A$ can be expanded as $A=\sum_i \lambda_i \Pi_i$, where $\lambda_i$ is the $i^{th}$ eigenvalue and $\Pi_i$ the orthogonal projection operator onto the $i^{th}$ eigenspace. With this in mind, a ray $\boldsymbol \Psi$ in the underlying Hilbert space defines a probability distribution as follows. Pick any vector $\psi$ from that ray; then the probability of measuring $A$ to have value $\lambda_i$ is simply
$$\frac{\langle \psi, \Pi_i \psi\rangle}{\langle\psi,\psi\rangle}$$
With this as motivation, we define an event as an orthogonal projector $\Pi$, and the probability of that event to be
$$P_\Psi(\Pi):= \frac{\langle \psi,\Pi \psi\rangle}{\langle \psi,\psi\rangle} \in [0,1]$$
In that sense, $P_\Psi$ constitutes a probability measure on the set $\mathcal P(\mathcal H)$ of orthogonal projection operators. We now ask whether it is the most general form of probability measure; this question is answered in the negative by Gleason's theorem, which says that the set of all probability measures is in one-to-one correspondence with the set of positive-semidefinite, self-adjoint operators $\rho$ with $\mathrm{Tr}(\rho)=1$ (i.e. density operators), with the corresponding probability measure being given by $\Pi \mapsto \mathrm{Tr}(\Pi \rho)$. The probability measures which arise from rays in $\mathcal H$ - which correspond to so-called pure states - form a strict subset of all probability measures.
The proof of Gleason's theorem is technical and beyond the scope of a PSE answer in my opinion. However, it is useful to see how the set of pure states is insufficient. Consider $\mathcal H= \mathbb C^2$ - does there exist a state of the system whose probability measure is uniform, such that for any observable $A$ with two distinct outcomes (eigenvalues) the probability of each is 1/2?
If we restrict ourselves to pure states, then the answer is no which can be seen easily. Consider $P_\Psi$ for some arbitrary ray $\Psi$, let $\psi$ be a normalized vector in that ray, and let $\phi$ be a normalized vector in the orthogonal complement of $\Psi$. Then $\{|\psi\rangle\langle\psi|,|\phi\rangle\langle\phi|\}$ is a set of orthogonal projectors which sum to the identity, yet $P_\Psi(|\psi\rangle\langle\psi|)=1\neq 1/2$ and $P_\Psi(|\phi\rangle\langle\phi|) = 0\neq 1/2$.
As a result, if we want to generalize to allow for such statistical ensembles, we must enlarge the state space from the set of rays in $\mathcal H$ to the set of density operators on $\mathcal H$. It's not hard to show that the latter objects can be written $\rho = \sum_i p_i |\psi_i\rangle\langle \psi_i|$ where $\sum_i p_i = 1$ and $\{\psi_i\}$ is some collection of unit vectors in $\mathcal H$.