Fatigue lives of rolling element bearings exhibit a wide scatter due to the statistical nature of the mechanisms responsible for the bearing failure process. Life models that account for this dispersion are empirical in nature and do not provide insights into the physical mechanisms that lead to this scatter. One of the primary reasons for dispersion in lives is the inhomogeneous nature of the bearing material. Here, a new approach based on a discrete material representation is presented that simulates this inherent material randomness. In this investigation, two levels of randomness are considered: (1) the topological randomness due to geometric variability in the material microstructure and (2) the material property randomness due to nonuniform distribution of properties throughout the material. The effect of these variations on the subsurface stress field in Hertzian line contacts is studied. Fatigue life is formulated as a function of a critical stress quantity and its corresponding depth, following a similar approach to the Lundberg–Palmgren theory. However, instead of explicitly assuming a Weibull distribution of fatigue lives, the life distribution is obtained as an outcome of numerical simulations. A new critical stress quantity is introduced that considers shear stress acting along internal material planes of weakness. It is found that there is a scatter in the magnitude as well as depth of occurrence of this critical stress quantity, which leads to a scatter in computed fatigue lives. Further, the range of depths within which the critical stress quantity occurs is found to be consistent with experimental observations of fatigue cracks. The life distributions obtained from the numerical simulations are found to follow a two-parameter Weibull distribution closely. The life and the Weibull slope decrease when a nonuniform distribution of elastic modulus is assumed throughout the material. The introduction of internal flaws in the material significantly reduces the life and the Weibull slope. However, it is found that the Weibull slope reaches a limiting value beyond a certain concentration of flaws. This limiting value is close to that predicted by the Lundberg–Palmgren theory. Weibull slopes obtained through the numerical simulations range from 1.29 to 3.36 and are within experimentally observed values for bearing steels.