Misconception about the expectation of a quantum system For a two-level quantum system with energy eigenstates $|\phi_1\rangle$ and $|\phi_2\rangle$ at finite temperature, we can write a general state as $$|\Phi\rangle=c_1|\phi_1\rangle+c_2|\phi_2\rangle.$$
In quantum mechanics, we know that the expectation of an operator $\hat{O}$ is given by  $$\langle\Phi|\hat{O}|\Phi\rangle = |c_1|^2\langle\phi_1|\hat{O}|\phi_1\rangle +|c_2|^2\langle\phi_2|\hat{O}|\phi_2\rangle +c_1^*c_2\langle\phi_1|\hat{O}|\phi_2\rangle +c_2^*c_1\langle\phi_2|\hat{O}|\phi_1\rangle.$$
However, in statistical mechanics, the expectation value is instead
$$\text{probability}(\phi_1)\langle\phi_1|\hat{O}|\phi_1\rangle + \text{probability}(\phi_2)\langle\phi_2|\hat{O}|\phi_2\rangle,$$
which simplifies to only the first two terms of the expression obtained from quantum mechanics. What is wrong here?
 A: Your formula for the expectation value in statistical mechanics is wrong.
A pure quantum mechanical state in this system is a sum $$|\Phi\rangle = c_1 |\lambda_1\rangle + c_2|\lambda_2\rangle$$ of two basis vectors, with complex coefficients.   Note that I'm using your notation $\Phi$ for a general state, and switching to $|\lambda_1\rangle, |\lambda_2\rangle$ for the basis vectors.  If we have a list of states $|\Phi_i\rangle$, we get a list of coefficients $c_{i,1}$ and $c_{i,2}$.
A thermal state in this system is an expectation value coming from an arbitrary sum (or integral) of pure quantum mechanical states:
$$
E(\hat{O}) = \sum_i d_i\langle \Phi_i | \hat{O} | \Phi_i \rangle.
$$
Here the coefficients $d_i$ are real. 
You're getting confused because you're trying to equate a special two-state case of the latter formula with the expectation value coming from a pure state.  If you want to write the stat mech expectation value in terms of $|\lambda_1\rangle$ and $|\lambda_2\rangle$, the right formula is
$$
E(\hat{O}) = \sum_i d_i\big(|c_{i,1}|^2\langle\lambda_1|\hat{O}|\lambda_1\rangle +|c_{i,2}|^2\langle\lambda_2|\hat{O}|\lambda_2\rangle + c_{i,1}^*c_{i,2}\langle\lambda_1|\hat{O}|\lambda_2\rangle) +c_2^*c_1\langle\lambda_2|\hat{O}|\lambda_1\rangle\big)
$$
You get one kind of uncertainty from having a pure quantum state, and another from having a probability distribution on the space of quantum states.
A: Nothing is actually wrong.
The quantum answer includes $c_2^*c_1$ term which depends on the phase difference between $c_1$ and $c_2$.  But the phase difference is something that would vary randomly with even the tiniest coupling between states 1 and 2.  And thermodynamics requires such a coupling to exist to be applied.  So the stat mech answer "averages over" things that very randomly in thermal systems.  Since that phase varies, $$\langle c_2^*c_1\rangle = 0$$ averages out to zero, and the only terms that do not average out to zero in thermal equilibrium are the terms that both stat mech and quantum preserve.  
The quantum answer applies when the states have no interaction, in which case it is not correct to apply stat mech because stat mech gives the answer "after we have reached thermal equilibrium" and the smaller the interaction between 1 and 2 is, the longer it takes to reach equilibrium, and when the interaction is vanishingly small, the time to thermal equilibrium approaches infinity.  
A: Both of your formulae are correct. One critical point is that, in quantum mechanics, we are adding states with a particular phase. So, the combinations $|c_1\rangle + |c_2\rangle $ vs $|c_1\rangle + i|c_2\rangle $ vs $|c_1\rangle - |c_2\rangle $ are not terribly different, they're just different ways to combine states. In contrast, in statistical mechanics, you are adding real numbers (expectations values), so you had better use a $+$ and not a $-$.
It is also possible to do both at once. Take a look at the density matrix, which is a great way to handle that. There is a good introduction in Sakurai's Modern Quantum Mechanics, and this site looks decent.
Foregoing the density matrix formalism, this is how you would go about combining the two ideas. Take some states $|\Phi_i\rangle$ for $i=1,2$ (we'll stick to two levels) with a $[c_{ij}]$ matrix
$$ |\Phi_i\rangle = c_{i1} |\phi_1\rangle + c_{i2}|\phi_2\rangle  $$
This has an expectation value
$$ (A_i)_{quantum} = \langle \Phi_i | A | \Phi_i \rangle = (\sum_j c_{ij} \langle \phi_j |) A (\sum_j c_{ij} | \phi_j \rangle)$$
Now, we want to calculate some statistical quantity, and we say that our system has a weight $w_i$ of each quantum state, which has nothing to do with $c_{ij}$. For thermal samples, $w_i$ is proportional to the Boltzmann weight of the energy eigenstate $i$.
$$\langle A \rangle_{statistical} = \sum_{i} w_i \langle \Phi_i | A | \Phi_i \rangle = \sum_i w_i (A_i)_{quantum} $$
You might wonder how it would look in a different basis, that this, what weighting factors to use for $|\phi_i\rangle$. That is, which $v_i$ would I use to get the same answer for
$$\langle A \rangle_{statistical} = \sum_{i} v_i \langle \phi_i | A | \phi_i \rangle$$?
I'll leave it as an excercise for the reader. =)   Suffice to say, this transformation is one thing density matrices are useful for.
