Every singlet state of a quantum spin 1/2 system can be decomposed into a linear combination of valence bond basis states. The range of valence bonds within this linear combination as well as the correlations between them can reveal the nature of the singlet state, and are key ingredients in variational calculations. In this work, we study the bipartite valence bond distributions and their correlations within the ground state of the Heisenberg antiferromagnet on bipartite lattices. In terms of field theory, this problem can be mapped to correlation functions near a boundary. In dimension d >= 2, a non-linear sigma model analysis reveals that at long distances the probability distribution P(r) of valence bond lengths decays as |r|^(-d-1) and that valence bonds are uncorrelated. By a bosonization analysis, we also obtain P(r) proportional to |r|^(-d-1) in d=1 despite the different mechanism. On the other hand, we find that correlations between valence bonds are important even at large distances in d=1, in stark contrast to d >= 2. The analytical results are confirmed by high-precision quantum Monte Carlo simulations in d=1, 2, and 3. We develop a single-projection loop variant of the valence bond projection algorithm, which is well-designed to compute valence bond probabilities and for which we provide algorithmic details.