In this work we present a method for generating random matrices describing electromagnetic scattering from disordered media containing dielectric particles with prescribed single particle scattering characteristics. Resulting scattering matrices automatically satisfy the physical constraints of unitarity, reciprocity and time reversal, whilst also incorporating the polarization properties of electromagnetic waves and scattering anisotropy. Our technique therefore enables statistical study of a variety of polarization phenomena, including depolarization rates and polarization-dependent scattering by chiral particles. In this vein, we perform numerical simulations for media containing isotropic and chiral spherical particles of different sizes for thicknesses ranging from the single to multiple scattering regime and discuss our results, drawing comparisons to established theory.