When employing Neural Networks as trial functions for solving Partial Differential Equations, solutions are obtained by minimizing loss functions defined through integrals. Unlike trial spaces in Finite Element Methods (FEM), exact quadrature rules are unavailable, making integration errors unavoidable. In particular, Gauss-type rules may yield incorrect results analogous to overfitting. First, we analyze the impact of bias in quadrature choice. We demonstrate that using biased stochastic quadrature rules can yield incorrect solutions, even when they are asymptotically exact. To remedy this, we utilize unbiased rules. Second, we focus on variance reduction, which is essential for rapid convergence. We propose new unbiased Gauss-type rules for triangular and tetrahedral elements, giving greater flexibility in complex geometries. They are exact for polynomials of a prescribed order and substantially improve convergence rates compared to standard Monte Carlo methods. Finally, we discuss a quasi-interpolation property that yields automatic variance reduction as the neural network approaches the solution. We compare discretizations that satisfy this property (e.g., Physics-Informed Neural Networks and First-Order System of Least Squares) against those that do not (e.g., the Deep Ritz Method and Variational PINNs), observing that the former attain solutions that are an order of magnitude more accurate at the same computational cost.