How to prove that sum converges to integral using density of states? Essentially, I would like to prove
$$ \sum_k f(k) \to \int f(k) \rho dE \tag{1}$$
where
$$ \rho = \frac{dk}{dE} \tag{2}$$
is the density of states and $k \to \infty$.
The model is that there is a system with $k$ energy levels with energies $E_k$. We consider a limit with infinitely many energy levels ($k \to \infty$). We can assume that in this limit $f(k)$ becomes a continuous function.
I know that Riemann sums converge to the integral, for example
$$ \sum_i f(x_i) (x_{i+1} - x_i) \to \int f(x) dx \tag{3}$$
for sufficiently good function $f(x)$ and properly chosen partition (i.e., $\max(x_{i+1} - x_i) \to 0$). But I can't reduce my equation $(1)$ to Riemann integral definition.
 A: 
Essentially, I would like to prove
$$ \sum_k f(k) \to \int f(k) \rho dE \tag{1}$$
where
$$ \rho = \frac{dk}{dE} \tag{2}$$
is the density of states and $k \to \infty$.

As mentioned in the comments, you need to introduce a measure on the LRS to get the dimensions to work out. To put it another way, your $f(k)$ on the LHS can't be the same as your $f(k)$ on the RHS.
Probably what you are trying to say is that you have some function $f_k$ defined for a discrete set of $k$ (e.g., $k_1$, $k_2$, etc), and some other function $f(k)$ defined for a continuous variable that has the same values as $\frac{f_k}{\delta k}$ when evaluated at the discrete points, and where $\delta k = (k_{i+1}-k_i)$.
Then
$$
\sum_k f_k = \sum_k f(k)\delta k\;,
$$
and now we just say that f(k) is slowly varying enough such that we can pretend that $\delta k$ is small and say that the sum is approximately the same as the integral.
$$
\sum_i f_{k_i} = \sum_i f(k_i)\delta k \approx \int f(k) dk \;,
$$
I.e., we assume that
$$
\int_{\tt{within }\,\, \delta k\,\, \tt{of }\,\, k_i} dk f(k) \approx f(k_i)\delta k
$$
Then change variables to $E$ instead of $k$:
$$
\sum_i f_{k_i} = \sum_i f(k_i)\delta k \approx \int f(k) dk = \int f(k(E))\frac{dk}{dE}dE\;,
$$
