When I have was taught quantum statistical mechanics, the course derived the canonical ensemble by assuming that eigenstates of the total system Hamiltonian have equal a priori probabilities, and then imposing Lagrangian constraints for the average energy and sum of probabilities being equal to one.
I want to know, why physically it is the eigenstates of the Hamiltonian that should have equal a priori probabilities and not some other complete set of states that span the Hilbert space?
One way to derive the canonical ensemble is by maximizing the entropy of the subsystem subject to given constraints (usually including a constraint on the average total energy). In this derivation, no a priori probabilites are assumed at all. On the contrary, maximizing the entropy amounts to finding the least presumptuous probability distribution compatible with the given constraints.
Another way to derive the canonical ensemble is to consider a subsystem of a larger system that is in the microcanonical ensemble. The microcanonical ensemble consists of a subspace (of the full Hilbert space) specified by some set of constraints, which usually includes a constraint on the total energy $E$, such as $E<E_0$ for some given maximum $E_0$. The microcanonical ensemble assigns equal probabilities to all states in any orthogonal basis that spans the subspace. Here's the key: Within the given subspace, assuming equal probabilities to the vectors in one orthonormal basis is equivalent to assuming equal probabilities to the vectors in any other orthonormal basis. That's because when $p_n=1/N$ for all $n$, the density matrix $$ \rho=\sum_{n=1}^N p_n|n\rangle\langle n| $$ is proportional to the identity matrix within the $N$-dimensional subspace spanned by the orthonormal state-vectors $|n\rangle$, and this is true no matter which orthonormal basis is used. Using a basis in which the basis vectors $|n\rangle$ are eigenstates of the Hamiltonian is often convenient, but it's not necessary; all predictions depend only on the density matrix and so are independent of which orthonormal basis is used.
The same ansatz is made in classical stat mech to derive the expression for probability of a configuration in terms of the partition function: configurations of equal energy are equally likely to be found.
If you want to ask why the classical ansatz: dynamically, systems tend to go to lower energy states. A ball rolls down a hill, charges attract/repel. Therefore there is a correlation between how likely a configuration is and its energy. It's natural to posit that equal energy configurations are equally likely to be found.
