How to reflect quantum uncertainty in statistical mechanics? How does the superposition principle of the quantum mechanics reconcile with the statistical mechanics? Because in deriving canonical distribution, I don't see anything related to superposition principle, but canonical distribution can be used to describe quantum systems, for example, a single quantum oscillator in thermal contact with a heat reservoir.So how to derive canonical distribution not from the eigenstates of the Hamiltonian but from the superposition states?
More generally speaking, how to reflect quantum uncertainty in statistical mechanics?
 A: There are many cases that you can use statistical mechanics for various reasons. Regarding your question about the ensemble, it is possible to use them even in a single system.
Let's consider the most familiar example, the ideal gas, in which various physical properties have be derived. Now, let's divide the system with volume $V$ into $N$ sub-system with volume $V/N$. The $N$ should not be too large such that $V/N$ can be considered as microscopically large. Hence, all of the sub-system have the similar (mean) macroscopic properties as the originally system, but now you have $N$ realization of an ensemble instead of a single system! Even for interacting system, it is usually possible to divide them into multiple sub-system, if the division is larger than correlation length. Each sub-system can now be considered as a realization of the ensemble, and statistics can be performed for this ensemble. 
On the other hand, if the system has the ergodic property. Then the system at different time may be considered as a realization of the ensemble, with sufficiently long time separated. The time averaged properties will be the same as the ensemble averaged. As an example, the temperature may fluctuate over very short time period, but a slow response thermometer can report the average temperature.
A: We can use statistical mechanics whenever the phase space distribution of the system is independent of time. This means that the distribution has to be a function of the constants of motion of the system which are $E, \vec L, \vec p$. Generally we have control over $E$ and not the others so we generally deal with systems whose phase space distribution is purely a function if $E$. 
Now the way we go about applying this is by first working out the allowed energy level for a single particle. Here we consider the total Hamiltonian (potential) that a particle sees and then find out the allowed energy states. Note that this is not the same as taking a system of particles and choosing one of them and treating rest as heat bath. No no. We are solving for possible energy levels for a single particle. Once we have that, we will fill them up with the total number of particles such that the total energy is equal to the desired value. 
And there could be many ways of filling it such that we get the right total energy. These are all the possible microstates. This method works when the interaction between particles is negligible, because we ignore interactions when we solve for single particle states. 
A: The answer to your question "how to put quantum uncertainty into statistical mechanics" is: with the density matrix.
Full description is in textbooks etc,(Feynman's Lectures on Statistical Mechanics, and von Neumann's Mathematical Foundations of Quantum Mechanics are good) but briefly you have a standard ensemble (canonical, microcanonical, or whatever) of $N$ systems with system $i$ in state $|\phi_i>$ (not necessarily eigenstates).
Introduce the density operator, $\hat \rho={1 \over N} \sum_{i=1}^N |\phi_i> <\phi_i|$.
That's the statistical mechanics bit. Now the quantum. Introduce a basis set of states $|\psi_j>$ which are eigenstates of some operator you're interested in (you have a choice), and the density matrix  is defined using these states.
$\rho_{jk}=<\psi_j| \hat \rho|\psi_k>={1 \over N} \sum_{i=1}^N<\psi_j |\phi_i> <\phi_i|\psi_k>$. 
The diagonal elements $\rho_{kk}$ then give the probability that if you measure a random member of the ensemble it will turn our to be in the $k$ eigenstate, and that includes the statistical probability, from the ${ 1\over N } \sum$, and the quantum probability from the $|<\psi|\phi>|^2$.   It can also be used to give expectation values for operators.
The off-diagonal elements contain information about the nature of the uncertainty. A classic example is an ensemble of electrons which are half spin up and half spin down, compared to an ensemble which are all spin sideways (aligned along, say, the x axis). In both cases, both diagonal elements are just $1 \over 2$, whereas the off diagonals are zero in the first case but not the second. 
Incidentally, the whole 'collapse of the wave function'  business is equivalent to setting all off diagonal elements of the density matrix to zero.
