$\newcommand\myderiv[3]{\left(\frac{\partial #1}{\partial #2}\right)_{#3}}$ Let's talk a bit about fairness. We used the fairness to find the probabilities of being in the various eigenstates, by assuming that the “fairest” distribution would prevail. If you bring two separate systems together and allow them to equilibrate, then you would expect that the net fairness would either remain the same or would increase. This sounds a little like entropy in the second law, in that the net entropy of system plus surroundings can increase or stay the same, but cannot decrease. The maximum value of the fairness for a given system (which is the value it will have in equilibrium) is its entropy.
Solving for maximum fairness
Let's look the maximum value of the fairness (a.k.a. entropy), which is $$\mathcal{F} = -k \sum_i P_i\ln P_i$$ $$U = \sum_i P_i E_i $$ $$P_i = \frac{e^{-\beta E_i}}{Z}$$ $$\mathcal{F}_\text{max} = -k_B\sum_i P_i \ln P_i$$ $$= -k_B \sum_i \frac{e^{-\beta E_i}}{Z} \ln \left(\frac{e^{-\beta E_i}}{Z}\right)$$ $$= -k_B\sum_i \frac{e^{-\beta E_i}}{Z} \left( -\beta E_i - \ln Z\right)$$ $$= k_B\beta \sum_i \frac{E_i e^{-\beta E_i}}{Z} + k_B \sum_i \frac{\ln Z e^{-\beta E_i}}{Z}$$ $$\mathcal{F}_\text{max} = k_B\beta U + k_B\ln Z$$ At this point, we may want to solve for $U$ again, to get yet another relationship for $U$: $$U = \frac{\mathcal{F}_\text{max}}{k_B\beta} - \frac{1}{\beta}\ln Z \qquad\qquad (9)$$ We saw before that $\ln Z$ was extensive, so we can now conclude that $\beta$ is intensive. From which it is also clear that entropy is extensive (which we already knew). Since we believe that $$S = \mathcal{F}_\text{max}$$ let's see what else we can extract from Equation 9 for $U$. We also know that $$dU = TdS - pdV$$ $$T = \myderiv{U}{S}{V}$$ Since we have an equation for $U$ in terms of $S$, we just need to figure out how to hold $V$ constant, and we'll know what $T$ is!
What does it mean to hold $V$ constant? It hasn't shown up in any of our statistical equations? If we change the volume, we will change the energy eigenvalues, so if we hold $V$ constant (and in general, do no work) then the energy eigenvalues are fixed. So until we explicitly add states with different volumes, then $V$ is held constant in our calculations, and thus we should be able to evaluate the derivative of $U$ with respect to $S$ at fixed $V$ to find the temperature. $$U = \frac{S}{k_B\beta} - \frac{1}{\beta}\ln Z$$ $$T = \myderiv{U}{S}{V}$$ $$= \frac{1}{k\beta} - \frac{S}{k\beta^2}\myderiv{\beta}{S}{V} + \frac{\ln Z}{\beta^2}\myderiv{\beta}{S}{V} - \frac1{Z\beta}\myderiv{Z}{S}{V}$$ $$= \frac{1}{k\beta} - \frac{S}{k\beta^2}\myderiv{\beta}{S}{V} + \frac{\ln Z}{\beta^2}\myderiv{\beta}{S}{V} - \frac1{Z\beta}\myderiv{Z}{\beta}{V} \myderiv{\beta}{S}{V}$$ $$= \frac{1}{k\beta} + \frac1{\beta}\left(- \frac{S}{k\beta} + \frac{\ln Z}{\beta} - \frac1{Z}\myderiv{Z}{\beta}{V} \right) \myderiv{\beta}{S}{V}$$ $$= \frac{1}{k\beta} + \frac1{\beta}\left(- \frac{S}{k\beta} + \frac{\ln Z}{\beta} + U \right) \myderiv{\beta}{S}{V} \qquad\qquad (18)$$ $$= \frac1{k\beta}$$ $$\beta = \frac1{kT}$$ where in step 18, we used the equation for $U$, Equation 9, which saved us the tedium of evaluating $\myderiv{\beta}{S}{V}$. By inserting this definition for $\beta$, into Equation 9, we can see that $$U = TS - kT\ln Z$$ $$-kT\ln Z = U - TS$$ $$F = -k_BT \ln Z$$ So it turns out that the log of the partition function just about gives us the Helmholtz free energy! $\ddot\smile$ This is often a bit more useful than our expression for $U$, since we know the derivatives of $F$ with regard to $T$: $$dF = -SdT - pdV$$ So we could conveniently evaluate $S$ by taking a derivative of the Helmholtz free energy—which would take us in a bit of a circle, but would allow us to express $S$ directly in terms of the partition function $Z$. Once we have computed $F(T,V)$ using statistical mechanics—which really only requires that we evaluate the partition function—we can use ordinary thermodynamics to compute all other thermodynamic quantities!