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 Hibbs 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}.
$$