Many modern outer radiation belt models simulate the long‐time behavior of high‐energy electrons by solving a three‐dimensional Fokker‐Planck equation for the drift‐ and bounce‐averaged electron phase space density that includes radial, pitch‐angle and energy diffusion. Radial diffusion is an important process, often characterized by a deterministic diffusion coefficient. One widely‐used parameterization is based on the median of statistical ultra low frequency (ULF) wave power for a particular geomagnetic index Kp. We perform idealized numerical ensemble experiments on radial diffusion, introducing temporal and spatial variability to the diffusion coefficient through stochastic parameterization, constrained by statistical properties of its underlying observations. Our results demonstrate the sensitivity of radial diffusion over a long time period to the full distribution of the radial diffusion coefficient, highlighting that information is lost when only using median ULF wave power. When temporal variability is included, ensembles exhibit greater diffusion with more rapidly varying diffusion coefficients, larger variance of the diffusion coefficients and for distributions with heavier tails. When we introduce spatial variability, the variance in the set of all ensemble solutions increases with larger spatial scales of variability. Our results demonstrate that the variability of diffusion affects the temporal evolution of phase space density in the outer radiation belt. We discuss the need to identify important temporal and length scales to constrain variability in diffusion models. We suggest that the application of stochastic parameterization techniques in the diffusion equation may allow the inclusion of natural variability and uncertainty in modelling of wave‐particle interactions in the inner magnetosphere.