Here is is one way to justify this within quantum mechanics, that goes by the name of eigenstate thermalization hypothesis:
For any system we can expand states $|\psi\rangle$ in terms of energy eigenstates $|n\rangle$:
$$|\psi\rangle = \sum_m C_m |n\rangle $$
Then the expectation of an observable $O$ is given by:
$$\langle O \rangle (t) = \sum_{mn}C_m^* C_n O_{mn} e^{i(\omega_m-\omega_n)t}.$$
The time average of this is given by:
$$ \overline{O}=\sum_{mn} C_m^*C_n O_{mn} I_{mn}$$
$$ I_{mn}=\lim_{T\to \infty}\frac{1}{T}\int_0^T e^{i(\omega_m-\omega_n)}{\rm d}t $$
The above is fully general and just a matter of definitions. We now need to make assumptions: first, that there are no degeneracies $\omega_m\neq \omega_n$ and so $I_{mn}=\delta_{mn}$. This is quite remarkable: the time dependence of the system has completely gone away in the long time limit. What has happened is that no matter what carefully arranged phase relationships you might start your initial state in are, the fact they evolve at different times means that they eventually dephase. This then gives the time average as:
$$\overline{O} = \sum_m |C_m|^2 O_{mm}$$
Now the question this raises is that surely this should depend on the initial state, that is, on the $|C_m|$'s. However, it turns out that for many observables in the energy basis in a chaotic system, the matrix elements are given by:
$$O_{mn}=O(E_m)\delta_{mn}+\text{ Some small random matrix}$$
where $O(E_m)$ is the microcanonical expectation value (essentially by definition of the microcanonical ensemble). The fact that $O_{mn}$ looks almost diagonal means the average does not depend on the $|C_m|$'s and instead just gives:
$$\overline{O} = O(E) \sum_{m}|C_m|^2 = O(E) $$
ie the microcanonical ensemble average at the expectation value of the energy. The assumptions we had to make were that the off-diagonal (in the energy basis) matrix elements were small. This can be justified somewhat by random matrix theory, in which it is a result that if you replace your Hamiltonian with a random Hamiltonian then the corresponding random part of your observable scales like $D^{-1/2}$ in the dimension of your Hilbert space.
References
A good review of this is 'From Quantum Chaos and Eigenstate Thermalization to Statistical Mechanics and Thermodynamics' - ArXiV page here