We consider the nonlinear Schrödinger equation with dispersion modulated by a (formal) derivative of a $\alpha$-Hölder continuous time-dependent function. Due to the highly oscillatory nature of the problem classical numerical methods face severe order reduction in non-smooth regimes $\alpha < 1$. In this work, we develop a new randomized exponential integrator based on a stratified Monte Carlo approximation which allows us to average the high oscillations in the problem and obtain improved error bounds of order $\alpha + 1/2$. In addition, the new approach allows us to treat a far more general class of modulations than the available literature. Numerical results underline our theoretical findings and show the favorable error behavior of our new scheme compared to classical methods.