Correlation functions in the O(n) models below the critical temperature are considered. Based on Monte Carlo (MC) data, we confirm the fact stated earlier by Engels and Vogt, that the transverse two-plane correlation function of the O(4) model for lattice sizes about L=120 and small external fields h is very well described by a Gaussian approximation. However, we show that fits of not lower quality are provided by certain non-Gaussian approximation. We have also tested larger lattice sizes, up to L=512. The Fourier-transformed transverse and longitudinal two-point correlation functions have Goldstone mode singularities in the thermodynamic limit at k --> 0 and h=+0, i.e., G_perp(k) = a k^{-lambda_perp} and G_parallel(k) = b k^{-lambda_parallel}, respectively. Here a and b are the amplitudes, k is the magnitude of the wave vector. The exponents lambda_perp, lambda_parallel and the ratio b M^2/a^2, where M is the spontaneous magnetization, are universal according to the GFD (grouping of Feynman diagrams) approach. Here we find that the universality follows also from the standard (Gaussian) theory, yielding b M^2/a^2 = (n-1)/16. Our MC estimates of this ratio are 0.06 +/- 0.01 for n=2, 0.17 +/- 0.01 for n=4 and 0.498 +/- 0.010 for n=10. According to these and our earlier MC results, the asymptotic behavior and Goldstone mode singularities are not exactly described by the standard theory. This is expected from the GFD theory. We have found appropriate analytic approximations for G_perp(k) and G_parallel(k), well fitting the simulation data for small k. We have used them to test the Patashinski--Pokrovski relation and have found that it holds approximately.