Average entropy of quantum system in bipartite pure state for finite temperature [I got halfway through writing this when I found the paper that answers the question in (essentially) the affirmative.  I'll post it anyways in case anyone is interested.]
Background: If a random vector $\psi$ (pure quantum state) is drawn from a $NM$-dimensional vector space for finite $N$ according to the Haar measure, then the entanglement entropy 
$$S(\rho_N) = \mathrm{Tr}[\rho_N \mathrm{log} \rho_N], \qquad \rho_N = \mathrm{Tr}_M[\psi\psi^\dagger]$$
across a tensor decomposition into $N$ and $M$ dimensional vector spaces (subsystems) is highly likely to be almost the maximum
$$S_{\mathrm{max}} = \mathrm{log}_2(\mathrm{min}(N,M)) \,\, \mathrm{bits},$$
 for any such choice of decomposition.  More precisely, if we fix $\alpha=M/N$ and let $N\to \infty$, then the fraction of the Haar volume of states that have entanglement entropy more than an exponentially small (in $N$) amount away from the maximum is suppressed exponentially (in $N$).  In physics this was known as Page's conjecture (later proved), and for mathematicians it is a simple consequence of the concentration of measure phenomenon. 
For any given Hermitian operator $H$ on the vector space acting as a Hamiltonian, we can assign an expected value of the energy $\langle H \rangle_\psi = (\psi^\dagger, H \psi)$ to a given vector $\psi$.  The Haar measure, then, can be considered a Gibbs probability distribution $p_\psi \propto e^{-\beta\langle H \rangle_\psi}$ in the infinite temperature limit ($\beta \to 0$), i.e. all states are equally likely.  The Page conjecture is exactly what you'd expect if a subsystem of a global pure state is maximally mixed, given its size.
My question: is there a way to extend this statement about typical entanglement entropies to cases where the probability distribution is one for finite temperature?  More physically, can we prove that the typical entanglement entropy of a system in a pure global state is exponentially close to
$$S_\beta = \mathrm{Tr}[\rho^{(\beta)}_N \mathrm{log} \rho^{(\beta)}_N], \qquad \rho^{(\beta)}_N = \mathrm{Tr}_M[\psi\psi^\dagger]$$
with high probability, according to the Gibbs distribution $p_\psi \propto e^{-\beta\langle H \rangle_\psi}$?
 A: The answer by Popescu et al is basically "yes", and is so close to my question that I suspect I read this paper and forgot it. 


*

*Sandu Popescu, Anthony J. Short, and Andreas Winter,
"Entanglement and the foundations of statistical mechanics"
Nature Physics 2, 754 - 758 (2006) doi:10.1038/nphys444 [Free PDF].


The precise statement has been proved not with a Gibbs distribution on global pure states, but with a constraint (i.e. constant weight positive on some subspace, and 0 outside it; "constant" being determined by the Haar measure).  When the constraint is global energy, this turns out to be effectively very close to a Gibbs distribution because the density of states is exponential in the energy, so that the distribution of energy is highly peaked around the corresponding temperature.  Their statement is significantly stronger than mine in that they prove the reduced state of the system is actually very close to the canonical thermal state, not just that it has the same entropy.
The statement, eq. (2) in the article, is that for arbitrarily small $\epsilon$ there exist
$$\eta = \epsilon +\frac{1}{2}\sqrt{\frac{d_s}{d_E^{\mathrm{eff}}}} \quad , \quad \eta' = 4e^{-C d_R \epsilon^2}$$
such that
$$\frac{V[\{ \vert \phi\rangle \in \mathcal{H}_R \vert D(\rho_S(\phi),\Omega_S)\ge\eta\}]}{V[\{ \vert \phi\rangle \in \mathcal{H}_R\}]} \le \eta'.$$
Here, $D$ is the trace distance, $\mathcal{H}_R$ is the subspace of the global Hilbert space satisfying the constraint, $\rho_S$ is the reduced state of the system, $d_S$ and $d_R$ are the dimensions of the system and the constrained subspace, and 
$$d_E^{\mathrm{eff}} = \frac{1}{\mathrm{Tr} \Omega_E^2} \ge \frac{d_R}{d_S}$$
is the effective size of the environment.  In these expression, $\Omega_S = \langle \rho_S \rangle_{\mathcal{H}_R}$ and $\Omega_E = \langle \rho_E \rangle_{\mathcal{H}_R}$ are the Haar-average reduced states conditional on the constraint.
