We prove that the flow of the discrete nonlinear Schrödinger equation (DNLS) is the mean field limit of the quantum dynamics of the Bose–Hubbard model for N interacting particles. In particular, we show that the Wick symbol of the annihilation operators evolved in the Heisenberg picture converges, as N becomes large, to the solution of the DNLS. A quantitative $$L^p$$ L p -estimate, for any $$p \ge 1$$ p ≥ 1 , is obtained with a linear dependence on time due to a Gaussian measure on initial data coherent states.