Obtaining $\beta = \frac{1}{k_B T}$ from first principles derivation I'm nearly at the end of this derivation but totally stuck so I'd appreciate a nudge in the right direction
Consider a set of N identical but distinguishable particles in a system of energy E. These particles are to be placed in energy levels $E_i$ for $i = 1, 2 .. r$. Assume that we have $n_i$ particles in each energy level. The two constraints we impose are that $\sum_{i}^{r}n_i = N$ and $\sum_{i}^{r}E_i n_i = E$. 
The number of microstates in a given macrostate is given by
\begin{equation}
\Omega = \frac{N!}{\prod_{i}^r n_{i}!}
\end{equation}
We want to maximize this and for ease of notation, we work with $\ln\Omega$ and we use Stirling's approximation ($\ln x! = x\ln x - x$) to obtain
\begin{equation}
\ln\Omega = N\ln N - N - \sum_{i}^{r}n_i\ln n_i - n_i
\end{equation}
Maximizing this function subject to the constraints $\sum_{i}^{r}n_i = N$ and $\sum_{i}^{r}E_i n_i = E$ is a classic Lagrange multiplier problem. We represent the undetermined multipliers to be $\alpha$ and $\beta$ for the two constraints and obtain
\begin{align}
\frac{\partial\ln\Omega}{\partial n_i} &= \alpha\frac{\partial n_i}{\partial n_i} + \beta\frac{\partial E_i n_i}{\partial n_i} \\ \nonumber
\ln n_i &= \alpha + \beta E_i \\ \nonumber
\therefore n_i &= e^{\alpha}e^{\beta E_i}
\end{align}
Now, we use the first constraint equation to determine $\alpha$. We get
\begin{align}
\sum_i^r n_i &= N \\ \nonumber
\sum_i^r e^{\alpha}e^{\beta E_i} &= N \\ \nonumber
e^\alpha &= \frac{N}{\sum_i^re^{\beta E_i}} \\ \nonumber
e^\alpha &= \frac{N}{Z}
\end{align}
We have introduced the partition function, $Z=\sum_i^re^{\beta E_i}$ in the last line. Next, we have the second constraint equation that determines $\beta$
\begin{align}
\sum_i^r E_i n_i &= E \\ \nonumber
\frac{\sum_{i}^{r} E_i e^{\beta E_i}}{\sum_i^r e^{\beta E_i}} &= \frac{E}{N} \\ \nonumber
\end{align}
I'm assuming I should somehow connect $E$ with $T$ so let's say $E=Nk_B T$. Then we have 
\begin{align}
\frac{N}{Z}\frac{\partial Z}{\partial\beta} &= E \\
\frac{\partial\ln Z}{\partial\beta} &= k_B T
\end{align}
How do I get to $\beta = -\frac{1}{k_B T}$ here? Notice that this derivation requires an extra minus sign compared to the usual definition of $\beta$ and this should come out naturally too, shouldn't it?
 A: I shall try to explain my understanding of how this problem should be solved. First, I agree with your derivation of
$$
n_i = N e^{\beta E_i}/Z, \qquad (1)
$$
where
$$
Z = \sum_{j=1}^r e^{\beta E_j} \qquad (2)
$$
and $\beta$ is obtained from equation
$$
\frac{E}{N} = \frac1Z \sum_{i=1}^r E_i e^{\beta E_i}. \qquad (3)
$$
We can make a conclusion now, that $n_i$ and $\beta$ depend on $E$ and $N$.
Our goal is to show that we have Gibbs distribution for $n_i$ and $\beta$ is equal to $-1/k_BT$. How can we relate $\beta$ and $T$? The only way I know is to use the following thermodynamic relation:
$$
\frac1T = \left(\frac{\partial S}{\partial E}\right)_N, \qquad (4)
$$
To obtain $S$ as a function of $E$ and $N$, we shall use Boltzmann's formula
$$
S = k_B\log\Omega.
$$
Then
$$
S = S_0(N) -k_B\sum_{i=1}^r n_i\log(n_i/e). \qquad (5)
$$
here $n_i$ depend on $E$ and $N$. Variation of $S$ induced by variations of $E$ and $N$ can be expressed through variations of $n_i$:
$$
\delta S = S_0'(N)\sum_{i=1}^r \delta n_i - k_B\sum_{i=1}^r\delta n_i \log(n_i). \qquad (6)
$$
Here $\delta n_i$ are expressed in terms of $\delta E$ and $\delta N$ in a somewhat cumbersome way. Substitution of (1) into (6) gives
$$
\delta S = -k_B\beta\sum_{i=1}^r E_i\delta n_i + \mbox{"smthng"}\sum_{i=1}^r \delta n_i. \quad (8)
$$
Obviously, $\delta n_i$ satisfy equations:
$$
\sum_{i=1}^r \delta n_i = \delta N, \qquad \sum_{i=1} E_i \delta n_i = \delta E.\quad (9)
$$
Then (8) and (9) give
$$
\delta S = -k_B\beta\delta E + \mbox{"smthng"}\delta N. \quad (10)
$$
The last equation gives
$$
\left(\frac{\partial S}{\partial E}\right)_N = -k_B\beta. \quad (11)
$$
At last, (11) and (4) give:
$$
\beta = -\frac1{k_BT}.
$$
Update. My explanation just demonstrates the statement of @higgsss's comment.
A: First of all, I think the number of microstates should be
$$\Omega=\frac{N!}{\Pi_in_i!}$$
$$\ln\Omega=\ln N! -\ln\Pi_in_i! $$
We have
$$\ln(\Pi_i x_i)=\sum_i\ln x_i $$
So
$$\ln\Omega=N\ln N -\sum_i\ln n_i! -N=N\ln N -\sum_i n_i\ln n_i+\sum_i n_i - N$$
$$\ln\Omega=N\ln N-\sum_i n_i\ln n_i$$
Otherwise I have no idea how did you simplify this:
$$\ln (\sum n_i!)$$
Or what is that supposed to mean. Check this out as to why I write $\Omega$ like that.
Now back to the problem. I can write
$$\sum_i n_i - N=0~~~and~~~\sum_i E_i n_i-U=0$$
$$\sum_i dn_i =0~~~and~~~\sum_i E_i dn_i=0~~~(*)$$ 
Let's take a look at differential of $\ln \Omega$
$$d\ln \Omega = d(N\ln N) -\sum \ln n_i dn_i -\sum n_i\frac{dn_i}{n_i}=-\sum \ln n_i dn_i$$
$$d\ln \Omega =-\sum \ln n_i dn_i=0$$
Now we can add (*) equations to the last equation, with two unknowns coefficients $\ln A$ and $\beta$:
$$-\sum \ln n_i dn_i-\ln A \sum_i dn_i +\beta \sum_i E_i dn_i =0$$
So
$$- \ln n_i -\ln A +\beta  E_i  =0$$
That's
$$n_i=e^A e^{-\beta E_i}$$
Let's find $A$ and $\beta$ to get it over with:
$$\sum n_i=N=\sum e^A e^{-\beta E_i}$$ 
And if 
$$Z=\sum e^{-\beta E_i}$$
$$e^A=\frac{N}{Z}$$
Next we should find $\beta$. Let's use this very important equation
$$\frac{1}{T}=\frac{\delta S}{\delta U}$$
First note that
$$\ln n_i  =+\beta  E_i -\ln A$$
So
$$d\ln \Omega =-\sum \ln n_i dn_i= -\sum (+\beta  E_i -\ln A) dn_i=-\beta \sum E_i dn_i -\ln A  \sum dn_i$$
The second term is obviously zero, but because $dE_i=0$ we can rewrite this equationt to get:
$$d\ln \Omega= -\beta d\sum E_i n_i=-\beta dU$$
right? (Of course $dU$ is zero in magnitude, but let's pretend we don't see that)
So if we consider the equation of entropy we have:
$$S=k_B \ln \Omega$$
And FINALLY
$$\beta = -\frac{d\ln \Omega}{dU}=-\frac{\delta S}{k_b\delta U}=-\frac{1}{k_bT}$$
My fingers hurt. Hope this help.
