Here is an explanation that's purely quantum.
A charged quantum particle in a magnetic field is subject to Landau quantization. Taking the magnetic field in the $z$ direction, we can choose the Landau gauge for the vector potential:
$$ \mathbf{A} = B x \hat{y} ~~ \Rightarrow ~~ \mathbf{B} = B \hat{z}. $$
The Hamiltonian in the coordinates $xy$, ignoring (for now) the edges of the sample:
$$ H = \frac{1}{2m} \left( \mathbf{p} - \frac{e \mathbf{A}}{c}\right)^2 = \frac{1}{2m} \left[ p_x^2 + \left(p_y- m \omega_c x\right)^2\right],$$
where $\omega_c = eB/mc$ is the cyclotron frequency.
After separation of variables we get the wavefunctions:
$$ \psi(x,y) = f_n ( x- k_y / m \omega_c ) e^{i k_y y},$$
where $f_n$ are the eigenfunctions of the simple harmonic oscillator ($n=0,1,2...$). The expectation values of $p_y$ and $x$ for this wavefunction are $\langle p_y \rangle =k_y$ and $\langle x \rangle =k_y / m \omega_c$, and the current along the $y$ direction is proportional to the generalized momentum in that direction:
$$ \langle I_y \rangle = \frac{-e}{m} \langle p_y - m \omega_c x\rangle = \frac{-e}{m} (k_y - m \omega_c \frac{k_y}{ m \omega_c} )=0.$$
As expected, we get zero current in the bulk of the sample.
Now let's imagine we are near the edge of sample on the negative side of the $x$ axis. This means the particle will feel a confining potential $U(x)$ that looks roughly like:
This potential will deform the wavefunction $f_n$ to a wavefunction that has more weight in the positive direction of $x$ than before, and then we'll get $\langle x \rangle > k_y / m \omega_c$, leading to:
$$ \langle I_y \rangle > 0, $$
i.e. edge current in the positive $y$ direction. Notice that this is the same direction predicted classically.
