Is there a thermal ensemble of pure quantum states? Thermal expectation values in a quantum system are given by
$$
\langle \mathcal O \rangle = \mathrm{Tr}\;e^{-\beta H} \mathcal O / \mathrm{Tr}\; e^{-\beta H}
$$
with $\beta$ the inverse temperature and $H$ the Hamiltonian. To see that this is a generalization of the classical formula, look in a basis where $H$ is diagonal.
But, there's another way we could conceivably generalize classical statistical mechanics: define a probability distribution $p(\psi)$ on the space of pure quantum states, and take expectation values with respect to it:
$$
\langle \mathcal O \rangle = \int d \psi\; p(|\psi\rangle) \langle \psi| \mathcal O | \psi\rangle
$$
(Clarification: by "the space of pure quantum states", I mean the full set of normalized vectors, without any requirement of linear independence.)
Does such a distribution $p(\cdot)$ exist? Well, yes, of course one does:
$$
p(|\psi\rangle) = \sum_n e^{-\beta E_n} \delta(|\psi\rangle - |n\rangle)
$$
Here the sum goes over the eigenstates of $H$ (forgive my notation in the delta function...). This is just the distribution that reproduces the first formula above. Boring!
Counting arguments suggest there should be many more. For a two-level system, the space of density matrices is three dimensional (or if I made a mistake counting, at most finite), whereas the space of possible $p(\psi)$ is infinite-dimensional.
So: what other natural distributions $p(\cdot)$ are there that give the correct thermal ensemble?
 A: What you call "another way we could conceivably generalize classical statistical mechanics" seems to be the direct application of the previous expression when $e^{-\beta H}$ or $\mathcal{O}$ are diagonal:
\begin{align*}\label{eq:pareto mle2}
\langle\mathcal{O}\rangle &=   \frac{\operatorname{Tr} e^{-\beta H} \mathcal{O} }{\operatorname{Tr} e^{-\beta H}}\\
&=\frac{\int d\psi d\phi \left<\psi |e^{-\beta H}|\phi\rangle \langle \phi|\mathcal{O}|\psi \right> }{\int d\psi \left<\psi |e^{-\beta H} |\psi \right> },\\
\end{align*}
which, under the assumption that $e^{-\beta H}$ is diagonal in the chosen basis, yields
$$\begin{align*}\label{eq:pareto mle23}
\langle\mathcal{O}\rangle &=\frac{\int d\psi \left<\psi |e^{-\beta H}|\psi\rangle \langle \psi|\mathcal{O}|\psi \right> }{\int d\psi \left<\psi |e^{-\beta H} |\psi \right> } \\
\end{align*}, $$
such that, defining
$$p(|\psi\rangle) = \frac{\langle \psi |e^{-\beta H}|\psi\rangle }{\int d\psi \left<\psi |e^{-\beta H} |\psi  \right>} $$
yields your result (you are missing a normalization factor). This is hardly an alternative definition.
Moreover, for a two level system there are only two states, implying that the density matrix is 4 dimensional. There are a finite number of states, hence $p(|\psi\rangle)$ is "two-dimensional" (although you can have an infinite number
of linear combinations of states).
What you asked, might be possible if the integral is extended beyond a basis of the Hilbert space. In that case there will be redundancy that would allow you to make  choices for $p(|\psi\rangle)$. But even in this case, these choices would only be possible for a given operator $\mathcal{O}$. Overall, that is not very useful.
So, to answer your question, there are no other distributions $p(\cdot)$, since the one you found is the original and there is no freedom to choose, as you have no freedom to choose the energy of a state.
For the case of two-dimensional system you have that:$$p(|1\rangle)= e^{-\beta E_1}\quad p(|2\rangle)= e^{-\beta E_2}.$$ There aren't any other choices that will yield the same results for all possible operators (at least for finite dimensional basis). (As an exercise you can try to find one).
Your question is fundamentally if there are two distinct matrices such that the trace of their product with any other matrix is the same:
$$\operatorname{Tr}{A\mathcal{O}}/\operatorname{Tr}{A} = \operatorname{Tr}{B\mathcal{O}}/\operatorname{Tr}{B} $$
which implies that
$$\sum_{ij}(A/\operatorname{Tr}{A} -B/\operatorname{Tr}{B} )_{ij}\mathcal{O}_{ji} = 0 \forall \mathcal{O}.$$
By restricting $\mathcal{O}$ to the hermitian matrices you have $n^2$ parameters. The above equation is of the form
$$\sum_{n}x_n a_n = 0\forall a_n$$
which can only be satisfied when $x_n = 0$, implying that $(A/\operatorname{Tr}{A}-B/ \operatorname{Tr}{B})_{ij}=0$. We conclude that there is a single matrix $A$ (up to multiplication by a constant) such that
$$
\langle\mathcal{O}\rangle =   \frac{\operatorname{Tr} e^{-\beta H} \mathcal{O} }{\operatorname{Tr} e^{-\beta H}}=   \frac{\operatorname{Tr} A \mathcal{O} }{\operatorname{Tr} A}.
$$
As the distribution $p(|\psi\rangle)$ can be obtained from $e^{-\beta H}$ and this quantity is uniquely defined, there is a single function $p(|\psi\rangle)$ that will yield the right result (up to multiplication by a constant).
For infinite dimensional spaces, the above condition will become
$$\int dx f(x)g(x) = 0 \forall g(x),$$ whose only solution is probably f(x)=0 (but I am not very familiar with the properties of infinite dimensional spaces, hence there might be some nuances that I am missing).
