Question: Newton's third law states that the force object A exerts on object B is equal to the force object B exerts on object A (in the opposite direction). So if, for example, a 5 kg ball is dropped from a plane, the force of the ball is 50 (F = 5\*10), so if the ball is exerting a force of 50N on the air molecules around it (even if it is one air molecule), shouldn't those molecules exert 50N on the ball too? And so shouldn't that ball always be in terminal velocity?

Answer: 
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.
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.