A hybrid finite element–boundary integral–characteristic basis function method (FE-BI-CBFM) is proposed for an efficient simulation of electromagnetic scattering by random discrete particles. Specifically, the finite element method (FEM) is used to obtain the solution of the vector wave equation inside each particle and the boundary integral equation (BIE) using Green's functions is applied on the surfaces of all the particles as a global boundary condition. The coupling system of equations is solved by employing the characteristic basis function method (CBFM) based on the use of macro-basis functions constructed according to the Foldy–Lax multiple scattering equations. Due to the flexibility of FEM, the proposed hybrid technique can easily deal with the problems of multiple scattering by randomly distributed inhomogeneous particles that are often beyond the scope of traditional numerical methods. Some numerical examples are presented to demonstrate the validity and capability of the proposed method.