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).
