La respuesta a tu pregunta "cómo introducir la incertidumbre cuántica en la mecánica estadística" es: con la matriz de densidad.
La descripción completa se encuentra en los libros de texto, etc, (Conferencias de Feynman sobre Mecánica Estadística, y von Neumann Fundamentos Matemáticos de la Mecánica Cuántica son buenos), pero en pocas palabras usted tiene un conjunto estándar (canónica, microcanónica, o lo que sea) de $N$ sistemas con sistema $i$ en el estado $|\phi_i>$ (no necesariamente estados propios).
Introducir el operador de densidad , $\hat \rho={1 \over N} \sum_{i=1}^N |\phi_i> <\phi_i|$ .
Esa es la parte de la mecánica estadística. Ahora la cuántica. Introducir un conjunto base de estados $|\psi_j>$ que son eigenestados de algún operador que te interese (puedes elegir), y el matriz de densidad se define utilizando estos estados. $\rho_{jk}=<\psi_j| \hat \rho|\psi_k>={1 \over N} \sum_{i=1}^N<\psi_j |\phi_i> <\phi_i|\psi_k>$ .
Los elementos diagonales $\rho_{kk}$ la probabilidad de que si se mide un miembro aleatorio del conjunto, éste se encuentre en la posición $k$ eigenestado, y que incluye la probabilidad estadística, a partir de la ${ 1\over N } \sum$ y la probabilidad cuántica del $|<\psi|\phi>|^2$ . También se puede utilizar para dar valores de expectativa de los operadores.
Los elementos no diagonales contienen información sobre la naturaleza de la incertidumbre. Un ejemplo clásico es un conjunto de electrones con la mitad de espín hacia arriba y la otra mitad hacia abajo, comparado con un conjunto con espín lateral (alineado, por ejemplo, a lo largo del eje x). En ambos casos, los dos elementos diagonales son simplemente $1 \over 2$ mientras que las diagonales son cero en el primer caso, pero no en el segundo.
Por cierto, todo el asunto del "colapso de la función de onda" equivale a poner a cero todos los elementos no diagonales de la matriz de densidad.