Time-domain electromagnetic simulations involving polydispersive media, such as biological tissue, are not straight-forward. One popular approach for representing the polarization is the Cole-Cole (1936) model, a heuristic generalization of the standard Debye (1929) model which itself corresponds to a first order linear auxiliary ODE. The Cole-Cole model corresponds to a fractional order ODE which presents computational challenges. We examine an alternative approach based on using the Debye model, but with a probability distribution of relaxation times. Recent advances in the generalized Polynomial Chaos method make the idea of a random polarization attractive. We will describe this modeling approach and then introduce a discretization of the system coupled with Maxwell's equations. We will present a stability and dispersion analysis of the overall method and numerically demonstrate convergence rates.