Drag force is essential to dense flows, but accurate and robust drag model is still an open issue. Direct numerical simulations of a shallow and a deep bubbling bed of moving rigid spheres have been carried out in the present work by using an immersed boundary method, and big data are produced. It turns out that the drag force in fluidized beds is typically underestimated by traditional drag models which depend only on the particle Reynolds number and the void fraction. With two additional parameters representing the velocity fluctuation and position fluctuation of particles introduced, a drag model based on the artificial neural network is developed. Given the complicated structure of this model, a simplified drag model is also formulated by directly fitting the samples. The drag force predicted by both models agrees excellently with the DNS data and is much more accurate than that predicted by existing models.