Some of the key motivators of quantum theory came from statistical mechanics. I know you have a broad range of stat mech backgrounds, so I will try to make this self-contained. I want to get to the point where you can understand the relevant mysteries, and how quantum mechanics solves them.
The philosophy of statistical mechanics is that real systems are complicated. At best you have an approximate description of them. For example, you model the hydrogen atom with the single-particle Schrodinger equation, but you did not include the fact that there are other atoms around, or photons. If your model was good, these things you missed are not going to mess anything up on the short term, but you will probably have trouble in the long term. Similarly, if you repeat your experiment over and over, you are likely to see some variation. All you can really be sure of is the things which are protected by conservation laws.
For example, if you have an isolated system, then the energy should be conserved. However, for a big system you typically don't know what state you are in, nor do you really care. Instead of talking about properties of a particular state, you talk about the ensemble properties of all the states with the same energy. This is known as the ``microcannonical ensemble". Remarkably, on a coarse grained scale, most of the states with a given energy have the same properties, so ensemble predictions are actually quite good. [There are precise statements of this such as the ``eigenstate thermalization hypothesis".]
The next step is that you promote this approximation to dogma. You define ``equilibrium" to mean that you have an ensemble of systems, and that all states of the same energy are equally likely. Thus the probability of any one state of energy $E$ is $p=1/\Omega\equiv e^{-S/k_B}$, where $\Omega$ is the number of state with energy $E$, and $S=k_B \log(\Omega)$. This function $S$ is known as the entropy. The constant $k_B$ is historical -- it plays no role except to change units.
All states with fixed energy $E$ are equally likely.
The next thing to consider is a composite, where our ``system" is in thermal contact with a ``bath". The system+bath is isolated, so we expect that any state of the composite system with energy $E=E_S+E_B$ will be equally likely. If the system and the bath are only weakly coupled, then the energy eigenstates are product states: $\Psi_T=\psi_S \phi_B$. The probablility of any given system state $\psi_i$ (with energy $\epsilon_i$) will be proportional to the number of bath states with energy $E-\epsilon_i$:
\begin{equation}
p_i\propto \Omega_{\rm bath}(E-\epsilon_i)=e^{S_{\rm bath}(E-\epsilon_i)/k_B}.
\end{equation}
We now assume that the bath is very big. Then we will find $E\gg \epsilon_i$, and we can Taylor expand the exponent
\begin{equation}
p_i\propto \exp\left(S_{\rm bath}(E)-\frac{1}{k_B}\frac{\partial S_{\rm bath}}{\partial E} \epsilon_i\right).
\end{equation}
The first part, $S_{\rm bath}(E)$ is the same for all system states $i$, so it just changes the proportionality constant. The only important property of the bath is the rate of change of entropy with respect to energy. We call this the inverse temperature of the bath
\begin{equation}
\beta=\frac{1}{k_B T}\equiv \frac{1}{k_B}\frac{\partial S_{\rm bath}}{\partial E}.
\end{equation}
If this was a statistical mechanics course I would then show that this definition of temperature has all the right properties and reproduces all of thermodynamics. This is a quantum mechanics course, so I will skip that, and just reiterate the important result
\begin{equation}
p_i\propto e^{-\beta \epsilon_i}.
\end{equation}
The proportionality constant is traditionally given as $1/Z$. From the fact that $\sum_i p_i=1$, we find
\begin{eqnarray}\label{can}
p_i&=&\frac{e^{-\beta \epsilon_i}}{Z}\\
Z&=& \sum_i e^{-\beta \epsilon_i}.
\end{eqnarray}
Although $Z$ is just the proportionality constant, it turns out that if you can calculate $Z$ as a function of other parameters (such as temperature), then there are tricks which you can use to calculate almost all thermodynamics quantities. We call the ensemble described by Eq.~(\ref{can}) the canonical ensemble.
State $i$, with energy $\epsilon_i$, is found with probability $p_i=e^{-\beta \epsilon_i}/Z$ where $Z=\sum_i e^{-\beta \epsilon_i}$.
We can continue with this trick. Suppose there is some other quantity which is conserved -- say particle number. The microcanonical ensemble will then involve $S(E,N)$. If we connect our system to a heat and particle bath, we find
\begin{eqnarray}\label{grandcan}
p_i&=&\frac{e^{-\beta \epsilon_i+\beta \mu n_i}}{\cal Z},
\end{eqnarray}
where
\begin{eqnarray}
\beta\mu&=&-\frac{1}{k_B}\frac{\partial S}{\partial N}\\
\beta&=&\frac{1}{k_B}\frac{\partial S}{\partial E}\\
{\cal Z}&=&\sum_i e^{-\beta \epsilon_i+\beta \mu n_i}.
\end{eqnarray}
The quantity $\mu$ is known as the chemical potential, and this is known as the grand canonical ensemble.
State $i$, with energy $\epsilon_i$, and particle number $n_i$ is found with probability $p_i=e^{-\beta (\epsilon_i-\mu n_i)}/{\cal Z}$ where
${\cal Z}=\sum_i e^{-\beta (\epsilon_i-\mu n_i)}$.
These expressions are the ``$F=ma$" of statistical mechanics. Everything else is just an application. I'll hand out a few exercises where you get to try this out.
Occupation Numbers
Lets take an example. Suppose we have a piece of silicon with an impurity in it. This might be a chip inside your cell phone. Suppose the impurity has a potential which can trap an electron. We want to know how likely it is to have an electron bound there.
The electrons in the chip provide a bath with some chemical potential $\mu$. The impurity state has some energy $\epsilon$. There are two states of our impurity, occupied, and unoccupied. The occupied state has energy $\epsilon$ and number of particles $1$, so it occurs with probability
\begin{equation}
p_{\rm occ}=e^{-\beta(\epsilon-\mu)}/{\cal Z}.
\end{equation}
The unoccupied state has zero energy and zero particles, so it occurs with probability
\begin{equation}
p_{\rm un}=\frac{1}/{\cal Z}.
\end{equation}
The grand partition function is chosen to make the sum of the probabilities equal $1$, so
\begin{equation}
{\cal Z}= 1+e^{-\beta(\epsilon-\mu)}.
\end{equation}
Here is a question for you: what is the average number of fermions in the state?
Suppose we instead had an attractive potential on the wall of a container of liquid helium. Imagine the impurity state again has some energy $\epsilon$. Since these are bosons you could put more than one of them on the site. If we neglect interactions, the energy of having $n$ atoms in the state is $n\epsilon$, giving a probability
\begin{equation}
p_n=e^{-n\beta(\epsilon-\mu)}/{\cal Z}.
\end{equation}
The partition function will be
\begin{equation}
{\cal Z}=\sum_n e^{-n\beta(\epsilon-\mu)}=\frac{1}{1-e^{-\beta(\epsilon-\mu)}}.
\end{equation}
Here is a second question for you: what is the average number of bosons in the state?
These two results are really important, as the results hold for any state -- not just an impurity. Of course, it neglects interaction, which can be a problem.