 Solution

Maximize the Gibbs entropy subject to normalization and the two mean constraints. Introduce multipliers $\eta ,A,B$ and differentiateStationarity gives ${\rho }_i\propto \exp [-(A{E}_i+B{N}_i)/k]$. Strict concavity of entropy gives the constrained maximum, on its support. Differentiating the maximum value with respect to the constrained means, or differentiating the entropy along the maximizing distribution and using ${\sum }_{i}d{\rho }_i=0$, gives $dS=AdE+BdN$. The thermodynamic identification is $A=1/T$, $B=-\mu /T$. Therefore