We derive equations for the strongly coupled system of light and dense atomic ensembles. The formalism includes an arbitrary internal-level structure for the atoms and is not restricted to weak excitation of atoms by light. In the low-light-intensity limit for atoms with a single electronic ground state, the full quantum field-theoretical representation of the model can be solved exactly by means of classical stochastic electrodynamics simulations for stationary atoms that represent cold atomic ensembles. Simulations for the optical response of atoms in a quantum degenerate regime require one to synthesize a stochastic ensemble of atomic positions that generates the corresponding quantum statistical position correlations between the atoms. In the case of multiple ground levels or at light intensities where saturation becomes important, the classical simulations require approximations that neglect quantum fluctuations between the levels. We show how the model is extended to incorporate corrections due to quantum fluctuations that result from virtual scattering processes. In the low-light-intensity limit, we illustrate the simulations in a system of atoms in a Mott-insulator state in a two-dimensional optical lattice, where recurrent scattering of light induces strong interatomic correlations. These correlations result in collective many-atom subradiant and superradiant states and a strong dependence of the response on the spatial confinement within the lattice sites