In this paper a methodology is presented to model uncertainties in the hysteresis law of ferromagnetic materials. The uncertainties may arise for example from manufacturing imperfections. A phenomenological type of model, known as the Duhem model, is introduced to model the hysteretic properties. The Duhem model is described via two material functions and one scalar parameter. In order to incorporate the uncertainties in the law the two material functions and the scalar parameter are modeled as a cross-correlated trivariate random field. The random field is discretized with a minimal number of random variables by the truncated Karhunen-Loève expansion. The presented stochastic Duhem model is used to compute the statistics of the hysteresis loss in a simple toroidal single phase transformer for illustration purposes.