The stopping site of the muon in a muon-spin relaxation experiment (mu+SR) is in general unknown. There are some experimental techniques that can be used to guess the muon stopping site, but they often rely on approximations and are not generally applicable to all cases. In this work, we propose a purely theoretical method to predict muon stopping sites in crystalline materials from first principles. The method is based on a combination of ab initio calculations, random structure searching and machine learning, and it has successfully predicted the MuT_ and Mu_BC stopping sites of muonium in Si, Diamond and Ge, as well as the muonium stopping site in LiF, without any recourse to experimental results. The method makes use of Soprano, a Python library developed to aid ab initio computational crystallography, that was publicly released and contains all the software tools necessary to reproduce our analysis.