Q:

Why should the equal a priori states be eigenstates of the Hamiltonian?

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?

A:

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.

