We give an efficient classical algorithm that recovers the distribution of a non-interacting fermion state over the computational basis. For a system of $n$ non-interacting fermions and $m$ modes, we show that $O(m^2 n^4 \log(m/\delta)/ \varepsilon^4)$ samples and $O(m^4 n^4 \log(m/\delta)/ \varepsilon^4)$ time are sufficient to learn the original distribution to total variation distance $\varepsilon$ with probability $1 - \delta$. Our algorithm empirically estimates the one- and two-mode correlations and uses them to reconstruct a succinct description of the entire distribution efficiently.