Thomson scattering (TS) of electromagnetic radiation by relativistic electrons is of great interest as an electron distribution diagnostics in high temperature plasma devices. A Monte Carlo approach to TS simulation has been implemented. It is based on first principles, i.e. on the integration of electron trajectories and on the computation and averaging of radiated Liénard–Wiechert fields. The treatment is relativistic. Nonlinear effects, arising from ultraintense lasers (both linear and circular polarizations having been considered) and/or from highly energetic electrons, are included from the outset. TS spectra from isotropic monoenergetic electron distribution functions and from the relativistic Maxwellian distribution have been extensively studied, both for low-intensity and for ultraintense incoming laser. A particularly simple quadratic dependence of S2(ω) on ω has been found in the case of monoenergetic distributions. Numerical results are checked against known analytical results for low-intensity lasers, and found to be in agreement with them. For high-intensity lasers, new effects on scattered spectra arise, like an overall red-shift or the appearance of higher harmonics. The numerical methods developed in this work are general enough to deal with anisotropic and/or non-Maxwellian distribution functions.