The positron, as the antiparticle of the electron, can form metastable states with atoms and molecules before its annihilation with an electron. Such metastable matter–positron complexes are stabilized by a variety of mechanisms, which can have both covalent and noncovalent character. Specifically, electron–positron binding often involves strong many-body correlation effects, posing a substantial challenge for quantum-chemical methods based on atomic orbitals. Here we propose an accurate, efficient, and transferable variational ansatz based on a combination of electron–positron geminal orbitals and a Jastrow factor that explicitly includes the electron–positron correlations in the field of the nuclei, which are optimized at the level of variational Monte Carlo (VMC). We apply this approach in combination with diffusion Monte Carlo (DMC) to calculate binding energies for a positron e+ and a positronium Ps (the pseudoatomic electron–positron pair), bound to a set of atomic systems (H–, Li+, Li, Li–, Be+, Be, B–, C–, O– and F–). For PsB, PsC, PsO, and PsF, our VMC and DMC total energies are lower than that from previous calculations; hence, we redefine the state of the art for these systems. To assess our approach for molecules, we study the potential-energy surfaces (PES) of two hydrogen anions H– mediated by a positron (e+H22–), for which we calculate accurate spectroscopic properties by using a dense interpolation of the PES. We demonstrate the reliability and transferability of our correlated wave functions for electron–positron interactions with respect to state-of-the-art calculations reported in the literature.