Obtaining $\beta = \frac{1}{k_B T}$ from first principles derivation

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.