Quantitative numerical analyses of interacting dilute Bose-Einstein condensates are most often based on semiclassical approximations. Since the complex-valued field-theoretic action of the Bose gas does not offer itself to standard Monte Carlo techniques, simulations beyond their scope almost exclusively rely on quantum-mechanical techniques, which are inherently limited to small particle numbers. Here we explore an alternative approach based on a Langevin-type sampling in an extended state space, known as complex Langevin (CL) algorithm. While the use of the CL technique has a long-standing history in high-energy physics, in particular in the simulation of QCD at finite baryon density, applications to ultracold atoms are still in their infancy. Here we examine the applicability of the CL approach for a one- and two-component, three-dimensional non-relativistic Bose gas in thermal equilibrium, below and above the Bose-Einstein phase transition. By comparison with analytic descriptions at the Gaussian level, including Bogoliubov and Hartree-Fock theory, we find that the method allows computing reliably and efficiently observables in the regime of experimentally accessible parameters. Close to the transition, quantum corrections lead to a shift of the critical temperature which we reproduce within the numerical range known in the literature. With this work, we aim to provide a first test of CL as a potential out-of-the-box tool for the simulation of experimentally realistic situations, including trapping geometries and multicomponent/-species models.