We consider ordinary differential equations (ODEs) with random parameters. We focus on Monte Carlo (MC) sampling for computing the statistics of some quantities of interest (QoIs) given by the solution of the ODE problems. We use the 4th order accurate Runge-Kutta (RK4) method as the deterministic ODE solver. We then develop a hybrid MC sampling method that combines RK4 with neural network models to efficiently compute the statistics of QoIs within a desired accuracy. We present several numerical examples to verify the accuracy and efficiency of the proposed hybrid method compared to classical MC sampling. The hybrid method that we develop can be applied to more complicated physical systems, for instance given by partial differential equations.