Why quantum mechanics uses density operators instead of probability distributions over state space? Whenever I try to get my head around mixed states I am referred to the notion of density operators. I think that density operators were introduced to represent mixed states as operators.
For what I read I see that a mixed state is just a bunch of pure states where each state has assigned some probability which all sums to one. Sometimes a mixed state is referred as ensemble in the sense of statistical ensemble.
Let us constrain to the finite dimensional case. The one we have when we deal with quantum computers. So we have the Hilbert space $H(\mathbb{C}^n$) and the pure states are modeled by rays in $\mathbb{C}^d$. Thus, the state space is modeled by projective space $\mathbb{C}P^{n-1}$.
Let us fixed an orthonormal base $|e_1\rangle, \dots |e_n\rangle$ in $\mathbb{C}^n$ and a pure state $|\phi\rangle$. If we do a measurement is this base, then according to Born rule the probability of outcome being $|e_j\rangle$ equals $|\langle e_j|\phi\rangle|^2$.
If we have a mixed state composed of $|\phi_1\rangle,\dots|\phi_k\rangle$ with probability of $|\phi_i\rangle$ being $p_i$ then we consider density operator $\rho=\sum_ip_i|\phi_i\rangle\langle\phi_i|$. And then if we want to compute probability of outcome being
$|e_j\rangle$, we (this I am not sure) compute $\langle e_j|\rho|e_j\rangle$ which gives us $\sum_ip_i |\langle e_j|\phi_i\rangle|^2$.
And the above is what bothers me. This whole idea of density operators feels for me like overkill.
The formula $\sum_ip_i |\langle e_j|\phi\rangle|^2$ looks just like marginal distribution (which is a special case of pushforwad measure). Thus, instead of density operators we should model mixed states as probability distributions over $\mathbb{C}P^{n-1}$ (if one want to be more formal one would consider $\sigma$-filed generated by borel sets on $\mathbb{C}P^{n-1}$). Given such probability distribution $P$, the probability that $|\phi\rangle$ is in $S\subset \mathbb{C}P^{n-1}$ is given by $\int_S dP(|\phi\rangle)$. In the case of discrete probability measure we obtain precisely statistical ensemble of pure states.
In such framework we could consider the probability of outcome being $|e_j\rangle$ as a function of pure states. I.e. $pm_j:\mathbb{C}P^{n-1}\to[0,1]$ given by formula $pm_j(|\phi\rangle)=|\langle e_j|\phi\rangle|^2$. And since $pm_j$ is measureble it may be treated as random variable and we can do with it all fancy probabilistic stuff. We can integrate it $\int_{\mathbb{C}P^{n-1}}|\langle e_j|\phi\rangle|^2dP(|\phi\rangle)$ and obtain something which is the direct generalization of what we obtained above (I mean $\sum_ip_i |\langle e_j|\phi_i\rangle|^2$).
We can push this framework even further. An arbitrary observable $A:\mathbb{C}^n\to\mathbb{C}^n$ as a linear operator gives rise (or rather factor to) mapping $\bar{A}:\mathbb{C}P^{n-1}\to\mathbb{C}P^{n-1}.$ And since $\bar{A}$ is measurable it pushes forward the measure. Thus, we can apply the above reasoning as well.
I believe that what I outlined generalizes to arbitrary Hilbert space (and to arbitrary quantum operations).
Is the above probabilistic framework consistent with the one defined in terms of density operators? If so, shouldn't whole quantum mechanics consider only pure states interesting?
 A: You can write density matrices as ensembles of pure states and ensembles
of pure states as density matrices. Most people view density matrices
as simpler since they are unique, while many different ensembles of states
can give the same density matrix. As a simple example of that, for a 2 state
spin-half problem, an ensemble of probability 3/4 with spin quantized
along $z$ and $1/4$ with spin quantized along $-z$ gives the same
density matrix and the result of all measurements will be the same
as an ensemble that has probability $1/4$ with spin quantized along $x$
$1/4$ along $-x$ and $1/2$ along $z$. The density matrix formulation
makes this clear where the ensemble average of pure states (either
of the two sets) does not,
but does give the correct results for all measurements.
Another advantage to density matrices is that it is easy to "integrate out''
degrees of freedom, i.e. to trace over them in the density matrix formalism
to obtain reduced density matrices for subsystems. You could, do the
equivalent calculation with an ensemble of pure states to get an ensemble
of pure states in the reduced Hilbert space, but this would essentially
require calculating the reduced density matrix using the standard
formalism, and then finding its eigenvalues and eigenvectors. The
eigenvalues would be the probabilities and the eigenvectors would be the
pure states in the reduced Hilbert space.
Addendum:
Here are the 2 cases I mentioned written out in detail.
I will use the usual up/down basis with the usual Pauli matrices.
The state quantized along $z$ is
$\left (\begin{array}{c} 1 \\ 0\\ \end{array}\right)$.
along $-z$
$\left (\begin{array}{c} 0 \\ 1\\ \end{array}\right)$.
along $x$,
$\frac{1}{\sqrt{2}}\left (\begin{array}{c} 1 \\ 1\\ \end{array}\right)$.
and along $-x$,
$\frac{1}{\sqrt{2}}\left (\begin{array}{c} 1 \\ -1\\ \end{array}\right)$.
The density matrices, $\rho_{s,s'} = \langle s|\psi\rangle\langle \psi|s'\rangle$, are
for quantized along $z$
\begin{equation}
\left (
\begin{array}{cc}
1 & 0\\
0 & 0\\
\end{array}
\right)
\end{equation}
along $-z$
\begin{equation}
\left (
\begin{array}{cc}
0 & 0\\
0 & 1\\
\end{array}
\right)
\end{equation}
along $x$
\begin{equation}
\frac{1}{2}
\left (
\begin{array}{cc}
1 & 1\\
1 & 1\\
\end{array}
\right)
\end{equation}
and along $-x$
\begin{equation}
\frac{1}{2}
\left (
\begin{array}{cc}
1 & -1\\
-1 & 1\\
\end{array}
\right) \,.
\end{equation}
Preparing a system with $\frac{3}{4}$ quantized along $z$ and $\frac{1}{4}$ along $-z$ gives
a density matrix
\begin{equation}
\frac{3}{4} 
\left (
\begin{array}{cc}
1 & 0\\
0 & 0\\
\end{array}
\right)
+\frac{1}{4}
\left (
\begin{array}{cc}
0 & 0\\
0 & 1\\
\end{array}
\right)
= \left (
\begin{array}{cc}
\frac{3}{4} & 0\\
0 & \frac{1}{4}\\
\end{array}
\right) \,.
\end{equation}
Preparing a system with probability $\frac{1}{4}$ along $x$, $\frac{1}{4}$ along $-x$ and
$\frac{1}{2}$ along $z$ gives a density matrix
\begin{equation}
\frac{1}{4}\cdot
\frac{1}{2}
\left (
\begin{array}{cc}
1 & 1\\
1 & 1\\
\end{array}
\right)
+
\frac{1}{4} \cdot
\frac{1}{2}
\left (
\begin{array}{cc}
1 & -1\\
-1 & 1\\
\end{array}
\right)
+\frac{1}{2}
\left (
\begin{array}{cc}
1 & 0\\
0 & 0\\
\end{array}
\right)
= \left (
\begin{array}{cc}
\frac{3}{4} & 0\\
0 & \frac{1}{4}\\
\end{array}
\right) \,.
\end{equation}
which is identical. All measurements on these two ensembles will give the same result.
There are infinitely many such ensembles. Any ensemble can be described by the eigenvectors
and eigenvalues of the density matrix. This example uses a density matrix that is
diagonal
in the usual $z$ basis.
